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Section 1 

Literature Review on Mechanical Reliability and Probabilistic Design 


Prof. Paul H. Wirsching 
University of Arizona 


February 1985 



1. INTRODUCTION 

A simple illustration of the basic problem of structural or mechanical 
reliability (or "probabilistic design") is provided in Fig. 1. The problem 
is to ensure, that the probability of failure of the cantilever beam is accept- 
ably small under the action of the stochastic load, Q(t). Assume that localized 
yielding defines failure. The probability of failure p^ is then the proba- 
bility of the event that the maximum stress, S, corresponding to the maximum 
load, Q, (assuming that dynamics is not important) exceeds the yield strength, 

R. 

P f = P(R - S) 

( 1 ) 

= P ( R - ^%) 
bii 

Assume all factors (R, Q, L, b, h) possess uncertainty and are modelled 
as random variables. It follows from the definition of the joint probability 
density function that 

D f = P(R - S) = f f x (x)dx (2) 

r 'u 

wnere, in oeneral, X is the vector of all design factors, and f v is the joint 
probability density function of the random design factors, f is the region 
of- failure, i.e., where R - S. 

Relative to the PSAM project, a similar formulation can be used to con- 
struct a cumulative distribution function (cdf) of a stress or a response, U. 
Assume that U is a function of several design factors. 


i i 

U 


= f (X) 




The cdf of U is defined as 

Py(u) = P(U - U) 


(3) 


(4) 
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Fig. 1 An illustration of a simple design problem in which all design factors 
(R, Q, L, b, and h) can be considered as random variables. 




By analogy, Fp(u) can be evaluated by Eq. 2 where Cl = (LI - u). 

In the general case, solution of the mul ti -dimensional integral of Eq. 2 is 
impossible "in practice. Development of "probabilistic design theory" (or 
"structural reliability") is directed towards practical solution to problems 
of the type of Eq. 1 for the purpose of (a) reliability assessment of existing 
designs, and (b) development of probability based design requirements. 

Presented here is a narrative sumnary of literature in mechanical relia- 
bility (probabi 1 istic design). Contributions are growing and this presentation 
' s not comprehensive. However, there is confidence that most of the impor- 
tant works are cited. The emphasis in this review is for application to the 
PS AM project. 

2. HISTORICAL NOTE 

The origin of modern probability theory dates back to the 17th century 
when an ardent gambler, Chevalier de Mere consulted the Franch mathematician 
Blaise Pascal (1623-1662) regarding a problem about a game of chance. Pascal, 
in turn, corresponded with Pierre Fermat (1601-1665). Subsequently, there was 
a rapid growth in interest in the mathematics of probability applied to comes 
of chance. Karl Gauss ( 1 777-1855) and L 'ierre Laplac< (1749-1827) were the 
i irst to find applications in other rields. 3ut serious interest in the syste- 
matic application of probabilistic and statistical methods to structural and 
mechanical design did not develop until the mid-1950's. 

A brief history of the development of the theory of structural reliability 
is presented in the text by Lind, Krenk, and Madsen [L5] and in the Ph.D dis- 
sertation of Kjerengtroen [K3]. Parts of the following are quoted from their 
work. The history of structural reliability goes back some 50-60 years. The 
first phase appears in retrospect as a very slow beginning. Early pioneering 
contributions included those of Forsell [F5] and Mayer [M2], and later Easier [Bl]. 
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M. Prot published several papers (in French) from 1936 to 1951 on Statistical 
distribution of stresses. And Weibull developed statistical theories of 
strength; his name is now associated with an extreme value distribution of 
minima [W1 , W2]. Later Puqsley [P6] and Johnson [Jl] gave comprehensive 
presentations on the theory of structural reliability and of economical desicn. 

The modern era of probabil istic mechanical design started after the 
Second World War. In October 1945, a paper entitled, "The Safety of Struc- 
tures" appeared in the proceedings of the American Society of Civil Engineers. 
This historical paper was written by A. M. Freudenthal, and the purpose of 
it was to "analyse the safety factor in engineering structures in order to 
establish a rational method of evaluating its magnitude." It was selected 
for inclusion with many discussions in the 1947 Transactions of the American 
Society of Civil Engineers [F6], The publication of this paper marked the 
genesis of structural reliability in the U.S. Most of the ingredients of 
structural reliability such as probability theory, statistics, structural 
analysis and design, quality control, existed prior to that time. Neverthe- 
less, Prof. Freudenthal was the first to put them together in a definitive 
and comprehensive manner. He continued, for many years, to be in the fore- 
front of structural reliability and risk analysis as well as fatigue and 
fracture studies. A sample of his significant publications on structural 
reliability and fatigue are provided in Refs. F7 and F8. Another landmark 
paper in structural reliability which began to formalize analysis was written 
by Freudenthal, Garrelts, and Shinozuka and was published in 1966 [F9T. 

During the 1960‘s there was rapid growth of academic interest in 
structural reliability theory. Classical theory became well developed 
and widely known through a few influential publications such as that of 
Freudentnal, Garrelts, and Shinozuka [F9], Pugsley [P5], Kececioglu and 
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Cormier [K2], Ferry-Borges and Castenheta [F2], and Haugen [H3]. However, 
professional acceptance was low for several reasons. Probabilistic design 
seemed cumbersome, the theory seemed intractible mathematically and numer- 
ically. Little data were available. Modelling error was unknown. And 
system structural safety analyses seemed extraordinari ly complex. 

Tne early 1 960 1 s were spent in the search to circumvent these difficulties. 

Turkstra [T2] presented structural design as a problem of decision making 
under uncertainty and risk. Lind, Turkstra, and Wright [L2] define the 
problem of rational design of a code as finding a set of best values of 
the load and resistance factors. Cornell [C2] suggested the use of a second 
moment format, and subsequently it was demonstrated that Cornell's safety 
index requirement could be used to derive a set of safety factors on loads 
and resistance. This approach related reliability analysis to practically 
accepted methods of design. It has been modified and employed in many 
structural standards. 

In the ensuing years seme serious difficulties with the second moment 
format were discovered in the development of practical examples. First, it 
was not obvious how to define a reliability index in cases of multiple random 
variables, e.g., when more than two loads were involved. More disturbingly, 
Ditlevsen [D2] and Lind [L3] independently discovered the problem of invariance. 
Cornell's index was not constant when certain simple problems were reformulated 
in a mechanical equivalent way. Several years were spent in the search of a 
way out of the dilemma without resolution. In the early 1970's, therefore, 
second moment reliability based structural design was becoming widely accepted 
although at the same time it seemed impossible to develop a logically firm 

basis for the rationale. 
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The logical impasse of the invariance problem was overcome in the early 
1 970 ' s . Hasofer and Lind [H2] defined a generalized safety index which was 
invariant to mechanical formulation. This landmark paper represented a turn- 
ing point in structural reliability theory. Contributions, proposed in recent 
years, are extentions of the Hasofer-Lind approach which are more sophisticated 
mathematically. The era of modern probabilistic design theory which extends 
from the early 1970's to the present is reviewed in Section 6. 

3. GENERAL REFERENCE TEXTS ON PROBABILITY THEORY AND MATHEMATICAL STATISTICS 

There are a large number of text books on the market with the approxi- 
mate title of, Probability and Statistics for Engineers and Scientists." 

Two which are recommended for easy reading and reference are those by 
Meyer [M3] and Hines and Montgomery [H5]. At a slightly higher level is 
the work of Bowker and Lieberman [B4], Texts on the same level (upper 
class undergraduate) written by engineers for engineering practice include 
those of Benjamin and Cornell [B3] and the two volumes by Ang and Tang [A2, A3], 
Intermediate texts on mathematical statistics include those by Freund [F1C], 

Mood and Graybill [M4], and Lindgren [L6]. These present advanced topics, 
but in a Torm which can be understood by engineers having some background. 

Designs are often selected on the basis of extreme loads stresses or 
strains. Extreme value theory is described in most of the above ret'erc-nces , 
e.g., both Mood and Graybill [M4] and Lindgren [L6] have discussions on order 
statistics. And and Tang's second volume has a cnapter on extreme value 
theory [A3]. The elementary text by Hahn and Shapiro [HI] on statistical 
models in engineering, has a good elementary description of extreme value 
theory. However, the most definitive work on this topic, although it is 
difficult reading, is the text by Gumbel published in 1958 [G4]. 
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4. MECHANICAL RELIABILITY AND PROBABILISTIC DESIGN TEXT BOOKS AND GENERAL 
REFERENCES 

There are a few text books which provide elementary information on 
basic probabilistic theory and application to mechanical design. These 
include texts by Kapur and Lamberson [Kl] and Siddall [S2]. Haugen has 
written two books on probabilistic mechanical design, the first was 
published in 1968 [H3] and the second appeared in 1980 [H4]. These 
texts provide a wealth of practical information, but all fail to 
provide comprehensive summaries of modern techniques of reliability 
analyses developed in the past ten years. 

Two other texts which are very useful for many engineering applications 
are those of Mann et al. [Ml] and Lipson and Sheth [L4]. The former focuses 
upon classical reliability models and is an excellent reference for general 
applications, ... but no design theory. Lipson and Sheth have much useful 

information not considered elsewhere. 

Advanced text books which treat modern reliability theory include that 
of Elishakoff [El], Leporati [Ll], and Ang and Tang [A3]. A new text, 
not yet published, by Lind, Krenk, and Madsen [L5] is an advanced work which 
summarizes the mathematical theory of structural reliability. At this time, 
however, perhaps the most highly regarded text is that of Thoft-Christensen 
and Baker [Tl]. 

Other references which provide general summaries of modern design 
theory include works by Shinozuka [SI], Wirsching [W5] in addition to 
CIRIA 63 report [R3], and the NBS report of Ellingwood et al. [E2], 
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5. CONFERENCE PROCEEDINGS AND PERIODICALS 

In recent years there have been a number of specialty conferences on 
structural reliability. The International Conference on Structural Safety 
and Reliability ICCOSAR is held every three years, but conference proceed- 
ings are not readily available. The 2nd International Conference on Code 
Formats in 1976 was a particularly productive one and its proceedings are 
somewhat of a classic [Dl]. The ASCE has sponsored a series of four specialty 
conferences since 1969 entitled, Probabilistic Mechanics and Structural 
Reliability. Proceedings are available through ASCE for the 1979 and 1984 
conferences [PI, P3], ASCE has also sponsored a specialty conference in 
1981 and has published conference proceedings entitled, Probabi 1 i stic 
Methods in Str uctural Engineering which contain some excellent summary 
articles [P2]. 

There is a new journal entitled, Structural Safety (pub! i shed by Elsevier) 
strictly dedicated to structural reliability. In addition, the civil engineering 
profession has been perhaps most active in the development of modern structural 
reliability concepts and the ASCE Journal of Engineering Mechanics and 
J ournal of Structural Engineering contain almost monthly articles on prob- 
abilistic design theory. Survey and theme articles published in the Journal 
of Structural Engineering include a literature review on structural safety 
published in 1972 [S4j, a series of six articles in 1974 [S5], a series of 
eight articles in load and resistance factor design in 1978 [G2], and a 
series of four articles in fatigue reliability in 1982 [FI]. Moreover, 
mSCE committees (e.g., the ASCE Administrative Committee on Structural 
Safety and Reliability, . . . and its five working committees) have spon- 
sored many technical sessions and produced numerous articles. 
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6. A NARRATIVE SUMMARY OF THE DEVELOPMENT AND IMPORTANT REFERENCES OF 
"MODERN" MECHANICAL RELIABILITY ANALYSIS. 


6.1 Mean Value First Order Second Moment Methods (MVFOSM ) 

The beginning of modern probabilistic design theory can be arbitrarily 
defined by the introduction of the safety index by Cornell in 1969 [C3] . 
First define the "failure function," Z, or "limit state," so that the event 
Z - 0 is failure. For the example of Eq. 1 


Z = R 


= r 6 QL 


(5) 


In general, Z will be a function of k random design factors, X^. An 
approximation to the mean value of Z, an< ^ standard deviation of Z, o^, 
can be derived using the first terms of a Taylor's series expansion. 


U7 


Z(u) 


a 


2 

I 



2 

r~, 

i 


( 6 ) 

(7) 


where u. and are the mean and standard deviations of respectively. 

■j is the vector of mean values, (as a rule of thumb, higher order terms are 
significant if the COV's of the variables are greater than 15 r =). 

The safety index defined by Cornell is 



( 8 ) 


In the special case where Z is linear in normal variates (or Z is a 
multiplicative function in only lognormal variates), the probability of 
failure is exactly 

P f » v(-5) (9) 
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For all other cases the probability of failure so computed is only 
approximate. In some cases, tne approximation is not bad; in other 
cases the error is very significant. 

In complicated problems it is often relatively easy to numerically 
estimate y z and a^. This may be particularly the case in the PSAM 
project in which the random variables are related to a computer algo- 
rithm. Using a simple straightforward perturbation technique, the 
derivatives of Eq. 7 can be computed. 

6.2 Normal ar.d Lognormal Formats 

Before consideration of advanced methods, the two special cases 
where £ defines the exact reliability should be noted. Assume that the 
failure function is of the form 


m 

Z = A + V a. X (10) 

0 i=l 1 1 

where all X. are normal and the A's are constant. Z is normal, and the 
-7 and - z of Eos. 6 and 7 are exact as is p f of Eq. 9. This is known as 
the "normal format." In many cases, Haugen has shown that mechanical 
design problems can be approximated by the normal format [H3, H4 ] . 

Another coninon form seen in design is the multiplicative function in 
which the failure event can be written as 

n v 

A !■ X 1 - 1 (11) 

i = l 1 

Here A and the a.'s are constant. By taking the log of both sides of 
Eq. 11, the linear form of Eq. 10 results, . . . along with the condition 
that failure is defined by Z - 0. And if all X^ are lognormal, then Eq. 9 
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given an exact form for p^. This "lognormal format" which forms the 
analytical foundation for LRFD [G2], is employed by Wirsching for fatigue 
reliability calculations [W5, W7, W8]. 

Because a closed form expression for probability results from the 
normal or lognormal format, assumptions should be made wherever practical 
in the PSAM project to produce this simplified form. 


f . 3 Advanced Reliability Methods; The Generalized Safety Index 

The failure function, Z, of Cornell is defined so that failure is 
the event that Z - 0. So defined, the algebraic formulation of Z is net 
unique. For example, it would be equally valid to write the Z of Eq. 5 as 


L 3 



( 12 ) 


A fatal flaw in the safety index of Eq. 8 is that i is not invariant to 

mechanical formulation of the failure function, e.a., the 2 of Eq. 8 would 

depend upon whether Eq. 5 or Eq. 12 was used. 

In 1973 Hasofer and Lind presented a new definition of the safety 

index data which overcame the lack of invariance problem [H2]. The scheme 
works like this. Each "basic design variable" is transformed by sub- 
tracting its mean, u^, and dividing by its standard deviation o.. The 
"reduced coordinate" x^ so defined has mean of zero and standard deviation 
of one. Upon substitution into the failure function, a new failure function 
Z.(x) is defined in terms of these reduced variables x. The Hasofer-Lind 

1 'v, -V 

(H-L) generalized safety index is defined as the minimum distance from the 
origin of the reduced coordinates to Z,(x) the failure function in these 

* a. 

reduced coordinates. So defined, the generalized H-L safety index gives 
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the same value of 5 as the case where the limit state is linear and the 
variables are normal. In other cases, the H-L index, e» will differ from 
the Cornell index. In summary, the H-L generalized safety index provides 
a measure of reliability which is invariant to the mechanical formulation 
of the failure function, . . . and gives the same value as the Cornell 
index special case of the normal format. 

The concept of the generalized safety index may play a key role in the 
PSAM project. An estimate of the probability of failure can be made by 
employing Eq. 9 above with the H-L index. Even though no distributional 
information is used in the H-L index, probability of failure so defined will 
provide a reasonable estimate to the actual probability of failure in many 
cases. More generally, Eq. 9 can be used to construct a distribution 
function of a response. 

Pelative to the PSAM program, it is important to note that estimates 
of the cumulative distribution function (cdf) of functions of random variables 
can be made employing 6. Let U be a function of several design factors 

U = f(X.) (13) 

The cdf of U is defined as 

F u (u) = P(U - u) (14) 

So that by analogy, the failure function is Z = U - u. To construct the 
cdf, several values of u must be chosen, and the computation for 6 repeated, 
but the computations are rapid by digital computer. Wu and Wirsching have 
demonstrated this technique on a low cycle fatigue problem [W6]. 
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6.4 Advanced Reliability Methods; Rackwi tz-Fiessler and Chen-Lind 

A principal limitation of the H-l approach is that distributional 
information, even if available, is not used in the analysis. In 1978 
Rackwi tz and Fiessler suggested a method which extends the Hasofer- 
Lind safety index concept to accomodate distributional information of 
the design factors [R 1 ] - Their method transforms non-normal distributions 
into "equivalent" normal distributions by adjusting the mean and standard 
deviation so that the distribution and density functions of the non-normal 
variables and the equivalent normal variables are equal at the design point. 
This scheme, an iterative algorithm which converges to a safety index, is 
described in references E2, L5, R3, SI, T1 , and W7. It has been demonstrated 
by Wu that probability of failure using the Rackwi tz-Fiessler 3 in Eq. 9 
produces, in many cases, surprisingly good estimates of the probability of 
failure lV.'9j • Using a digital computer, the calculations are very fast and 
efficient. Again, the R-F method can be employed to construct distribution 
functions of responses in complicated problems, and therefore may be very 
useful in certain aspects of the PSAM. All of these schemes are referred to 
as "fast probability integration methods" because they provide a very fast 
approximation io Eq. 2 . 

An extension of the R-F scheme was proposed by Chen and Lind [Cl], 

This method uses a three parameter equivalent normal distribution. It 
was anticipated that this method can produce more accurate methods of 

p f than does R-F, but Wu has shown this to be not always true [W9], In 
fact, Wu has developed another advanced method, which employs techniques 

of Rackwitz and Fiessler and Chen and Lind to produce a safety index, and 
extimates of probability of failure with significantly less error than 
R-F or C-L. 
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6.5 Comments on Methods Which Rely on Higher Order Approximations to the 
Limit State 

The procedures described above are referred to as "first order metnods" 

because the limit state is approximated as a straight line at the design 
point. But other advanced methods have been proposed. In general, relia- 
bility analysis can be performed by transforming the basic variables, X.. , 
to standard normal variables, as suggested by Rosenblatt [R5] 

x = ^‘ 1 [F(X)J (15) 

the inverse transformation is 

X = F' 1 [:(x)] (16) 

The inverse transformation is substituted into the original failure func- 
tion, Z, so that the transformed failure function, Z 1 , can be formulated 
and the safety index computed. Such a procedure is expected to produce 
more accurate values of p^, but the inverse transformation can be extremely 
compl icated. 

Improvements to the first order method, suggested by various authors, 
typically employ a higher order approximation of the limit state. For 
example, Citlevsen developed bounds by inscribing and circumscribing the 
limit state with rotational paraboloids [D4]. Horn and Price investiaated 
the error of the linear approximation by studying an approximating hyper- 
sphere with radius corresponding to the mean curvature at the design point [H7]. 
To avoid the arbitrariness of the choice of a suitable approximating limit 
state, Fiessler et al. investigated several possible forms of the quadratic 
limit state [ F 5 ] . Breitung derived an asymptotic formula for the probability 
of failure which considers curvatures in the limit state at the design po i- nt 
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[B5]. Tvedt derived two approximation formulas, both which model the failure 
surface with a parabolic surface at the design point, which also give an 
accurate estimate of p^ [T3]. 

Because a more accurate description of the limit state is used, these 
second order methods have the promise of consistantly producing better 
estimates of the probability of failure relative to first order methods. 
However, these schemes are are more complicated because the formulations 
are made on the transformed space and require second partial derivatives 
of the transformed limit state function. In the literature, only a few 
simple examples have been presented, e.g., Fiessler et al. [F5], Breitung 
[B5], Tvedt [T3]. The evidence is not entirely convincing that quadratic 

methods could produce consistantly accurate results relative to other methods. 

In summary, it is not obvious at this time that the much more compli- 
cated second order methods will be helpful in the PSAM project. Wu has 
shown that his first order method works extremely well, and errors in 
computing probabilities are almost always acceptably small [W9]. 

6.6 On the Drawing Board 

While it is widely recognized that higher order forms for approximating 
the limit state may produce more accurate estimates of probabilities, they 
are typically far more complex than linear forms. Wu has developed a linear 
limit state algorithm, an extension of the R-F and C-L methods, which is 
efficient and accurate [W9]. He also has under development at this time an 
advanced version which, based on a few check cases, seems to be faster and 
more accurate. 

6.7 Reliability Analysis When the Limit State Function Does Not Have a Closed 
Form Expression 

Reliability methods described above rely on a closed form expression of 

the limit state, e.g., Eq. 5. But there are cases where the relationship 
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between the variables is defined only by a computer algorithm, e.g., local 
strain and fracture mechanics fatigue life prediction, finite element stress 
analysis, etc. Wu and Wirsching have presented a method whereby the computer 
program is run using various combinations of parameter values, and a poly- 
nomial approximating the limit state is fitted to the responses [W6]. A 
fast probability integration method is then employed to estimate probabilities. 
The good news, ... no modification of the program which characterizes 
physical behavior is necessary. The bad news, . . . particularly with 
regard to the PSAM project, ... a separate analysis must be done on each 
response variable. 

7. APPLICATION OF PROBABILISTIC DESIGN THEORY TO DESIGN CODE DEVELOPEMNT 

Probably the largest effort in the U.S. to implement a reliability based 
design criteria was the LRFD (Load and Resistance Factor Design) program to 
revise the AISC specifications. Work on the development of the new AISC 
specifications started in 1969, and it was conducted by M. K. Ravindra and 
T. V. Galambos. A comprehensive summary of the theoretical development is 
presented in eight papers in the ASCE Journal of the Structural Division , 

Sept. 1973 [G2] ; historical summary is provided by Galambos [G3]. The 
proposed specifications [P4] are now open for public review and discussion. 
After final revision the new rules will be included as an alternate to the 
1978 AISC specifications. The lognormal format, employed in LRFD, may be 
useful for elementary reliability analyses in the PSAM project. In particular, 
a closed form expression for the probability distribution of response is 
possible when the random variables (all assumed to be lognormal) can be 
factored outside of the stiffness and mass matrices in the static linear 
case. Furthermore, the partial safety factor format of LRFD can be employed 
if a safety check expression is required. 
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Other well documented efforts to develop probability based design re- 
quirements are available. An excellent description of the review and revision 
of the National Building Code of Canada is provided by Siu , Parimi, and Lind 
[S3]. The NBS report by Ellingwood et al. recommends load factors and load 
.combinations compatible with loads in the proposed 198C version of American 
National Standard A58 [E2]. Load factors were developed using concepts of 
probabilistic limit states design. Both the Ellingwood report and the CIRIA 
63 report [R3] provide excellent and comprehensive summaries of techniques 
and applications of modern probabilistic design theory. Moreover, Ellingwood, 
et al. [E2] and Galambos and Ravindra [Gl] provide useful data summaries on 
material behavior (structural steel at room temperature). The purpose of the 
CIRIA 63 report was to review suitable methods for the determination of partial 
factors for use in limit state structural codes. Both reports detail the 
definition and process for computing the generalized safety index, a technique 
which may be very useful in all aspects of the PSAM project. Another very 
excellent reference is the Bulletin d ' Information 112, published by the Comite 
Europeen du Beton [R2]. Unfortunately, this document is not readily available, 
but it does contain major contributions from many of the pioneers of the 
development of probabilistic design theory. 

8. A NOTE ON MONTE CARLO METHODS 

Monte Carlo is employed very effectively to analyze complicated problems 
in probability theory, mathematical statistics, reliability, random process 
theory, etc. As a general rule, Monte Carlo analysis tends to be very costly 
relative to the accuracy of the results. Therefore, it is commonly used in 
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a research role to verify the performance of more efficient numerical methods. 
In the PSAM project it is not likely to be an effective design tool. 

Monte Carlo is particularly inefficient for mechanical reliability 
problems because accurate estimates of the small probabilities of failure 
require very large sample sizes. Efficiency can be improved by discrimina- 
tion in sampling or by extrapolating an empirical distribution function; 
but generally speaking, advanced reliability methods cited in Section 6 are 
far more efficient for the basic reliability problem. 

Monte Carlo seems more of an art than a science, and no complete work, 
for engineering application, seems to exist. Elementary concepts (how to 
sample from various distributions) are presented in Hahn and Shapiro [HI], 

Both Ang and Tang [A3] and Elishakoff [El] have chapters on Monte Carlo 
presenting engineering applications. Thousands of papers have been pub- 
lished which describe a wide variety of applications. For example, two 
elementary works, by this author, describe application to random process 
simulation for fatigue analysis [W4], and analysis of peak responses to non- 
stationary random forces [W3]. 
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ABSTRACT 


The notion of stochastic variables in structural analysis 
was introduced by the late Professor A.M. Freudenthal as early as 
in 1945. The goal has been to assess structural safety in a 
rational fashion. One cannot totally rely on the hypothetical 
deterministic assumptions with the pretention that the knowledge 
is complete and exact regarding material properties, geometry of 
components, and loading. Hence the ^Probabilistic ^Structural 
Analysis ^Method (PSAM) emerged in order to evaluate structural 
performance in real world situations. Along with the advent of 
digital computers the finite element method has established 
itself to be the singlemost versatile numerical tool for 
engineering calculation. Stochastic analysis on the response 
database furnished by a finite element scheme is then the most 
logical way to carry out relevant reliability calculations for 
engineers who are responsible to assure safe functionality of 
systems they analyze, design and construct. 

Quantitative estimation of failure apprehension can be 
obtained by considering stochasticity of both loading and 
structural description. The former aspect is treated in random 
vibration and will not be addressed here. Available finite 
element type formulations with random variables describing 
stiffness, mass and damping matrices due to uncertainities in 
boundary geometry, initial stress distributions, . material 
properties and bounday conditions, are reviewed in this report. 

Computational procedure for evaluating the design statistics 
(such as the means, variations, correlations, etc.) of mode 
shapes, resonant frequencies, buckling loads and non-linear 
dynamic respnses are summarized. A list of reference of 
important publications is furnished. Comments on outstanding 
issues and necessary research is also included herein. 
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1 . introduction 

Engineering systems are designed with a variety of materials 
and are shaped conveniently in order to perform certain 
functions. During its service life a system enounters many 
different static and dynamic loading conditions. The main 
concern that spans from a lay person to a competent designer is 
(a) whether the structure will survive, (b) how well the behavior 
of the structure would correspond to the required specifications 
and, (c) what are the chances of encountering undesirable 
circumstances such as cracks and excessive vibrations. Everyone 
is interested in the the overall rating of performance as well. 

We can immediately detect that the direction of these natural 
questions are both quantitative as well as qualitative in 
nature. If we consider the entire design procedure to be a 
decision making activity, then at each instance we are compelled 
to resolve a generic question. What is the chance that certain 
cr iterion will not be met during the life of the engineering 
system which is conceived on a design board? 

We immediately recognize that the problem in engineering 
design analysis is bifocal. First, we must recognize physical 
behaviors and secondly, we must examine the extent of our 
knowledge regarding these behaviors. In order to answer the 
first question, we axiomatize a mathematical model and quantify 
applicable physical laws. Then analysis is performed adhering as 
closely as possible to exact solutions. Unfortunately, even many 
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timple objects of engineering design analysis are so complex in 
geometry that continuum methods succeeded by analytic solutions 
reduces to nothing more than text book examples. Thus in 
practice, based on the knowledge of systems of rather simplified 
geometry, discrete (as opposed to continuum) methodologies are 
pursued where the solutions are arrived at in numerial steps 
(contrary to analytical methods with closed form expressions). 
Computational methods such as finite difference and finite 
element techniques thus emerged as very powerful numerical 
tools. With the advent of high speed digital computers, it 
became possible to carry out a large number of arithmetic 
operations leading to the success of those numerical methods 
appropriate for dynamic response computation as well as thermal 
analysis. Thus the partial differential equations of 
mathematical physics, which dictate the motion, thermal behavior, 
etc., are reduced to rather simplified solution of algrebaic 
equations. The finite element method, which is a means to 
spatially discretize the continuum operator that governs the 
field variables, became very popular since the material 
inhomogeneity, anisotropy, arbitrariness of boundary geometry 
could be easily incorporated in the numerical procedures. In 
essence, the answer to the first question can be summarized in 
terms of applicability of the conventional finite element method. 

However, the second question invokes a different branch of 
discipline altogether viz. probabilistic analysis and statistical 
computations. We have first hand experience that the design 
assumptions are quite empirical if not gross to some extent. In 
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reality we are dealing with partial, often quite incomplete and 
contaminated information regarding the structure and loading 
conditions. Hence it is quite legitimate to attempt to evaluate 
differences between the predicted and any possible realistic 
responses. Very naturally, concepts like mean values, standard 
deviation, probability distributions and exceedence (probability 
to exceed the allowable limits) arise within the selected 
numerical method i.e., the finite element method. Thus a 
conjugation of the finite element procedure (spatial 
descretization) with the probabilistic notion of analysis becomes 
ineviable in a rational design-analysis environment. 

In order to illustrate the aforementioned generalized (to 
some extent rather vague) discussion let us consider one of the 
most simple problems in structural mechanics. This will also 
facilitate the introduction of some definitions like random 
variables, stochiastic processes, etc. which are vital to the 
appreciation of the cited literature reviewed in the succeeding 
chapters . 



Fig. 1.1 Uniaxial Bar Problem 
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Consider in Fig. 1.1 a uniform bar of length L, depth D, 
width B, subjected to a constant axial stress a. The material 
will be taken to be homogenous thus the modulus of elasticity E 
will be considered to be constant. Suppose we are interested in 
the strain e and elongation U of the member. From basic stength 
of materials: 

e- - and U * eL * — (1.1) 

E E 

Now we shall ask a pertinent question regarding our 
confidence in the assumptions leading to the expressions of e and 
U in equation (1.1). The first set of questions will address the 
loading. How accurately do we know that a is uniform on the end 
surfaces? If there is a device which applies the force we can 
never be sure that a perfectly uniform stress condition is 
imposed. The rational way to proceed will be to estimate 
functions F ff (x,y , a) on the left and right faces such that at a 
point (x,y) the probability of the applied stress to be less than 
o will be given by the value of the function F. At this stage 
let us assume that the bar is “perfect" with its stipulated 
geometrical dimensions and modulus of elasticity. The resulting 
strain distribution e(x,y,z) will also now become uncertain as a 
consequence of the distribution F 0 (x,y,o). Then the pertinent 
design quantity to look for, in order to perform an analysis on 
the basis of strain, will be F e (x,y,z,c), i.e., the probability 
distribution function for the strain e. It is interesting to 
note that the stochastical strain now becomes a three-dimensional 
function even for the corresponding one-dimensional deterministic 
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case . 


Hence we need to carry out a three— diinensional analysis of 
the aforementioned bar of Fig. 1.1. In order to utilize an 
available finite element computer program we spatially discretize 
this static problem. Without any loss in generality and 
especially in order to avoid unnecessary complexity let us assume 
that the end stresses are so applied that a does not vary with 
x. Hence the probability distribution function for o could be 
represented in the form F e (y,o). If from our engineering insight 
we assume that the resulting strain e does not vary with x at a 
certain section then we would like to evaluate F £ (y.z,e). For 
this two-dimensional idealization we employ a two-dimensional 
finite element mesh as shown in Fig. 1.2. 



Fig. 1.2 Finite Element Discretization 
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The objective of the probabilistic finite element procedure within the 
context of the problem in Figs. 1.1 and 1.2 is to evaluate F E (y,z,e) 
when F (y,o) are given for each end faces. 

In the preceeding example we consider the stress, a, loading 
(forcing function in the finite element system) to be a 
nondeterministic function of x and y. This will be termed to be a 
stochastic or random process . A formal definition of a random process 
is that f(x) is a random process if f is a random function of a 
deteministic argument x. We shall indicate stochastical variables 
with a tilda. Thus a random function such as the stress, strain in 
the above problem will be random processes, c(y) and e(y,z), 
respectively. (Notation introduced in [B-1.4] will be used throughout 
this report . ) 

A dynamic version of the above class of problems, depicted in 
Figs. 1.1 and 1.2, attracted the notice of several researchers. 

Therein a structure or any other mechanical system was considered to 
be completely deterministic whereas only the forcing function (such as 
the earthquake or wind load) was considered to be random processes in 
time. This class of problem of deterministic system with random 
loading are treated in a special branch of dyamics called random 
vibration, refer to [C-1.4]. There are excellent standard text books 
on that topic as listed in the reference (section 10). 

We can further pursue our question regarding the assumptions in 
the formulation leading to t a | , U * in equation (1.1). There 

are possibilities that during manufacturing of the "real world" bar in 
Fig. 1 the chemical process was not exact or perfect hence the modulus 
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of elasticity E may just be almost constant and is indeed a random 
function of x,y.and z, i.e. E(x,y,z). Since the manufacturing 
processes are quite reliable we expect E to be a narrow band process 
implying that the difference in say the maximum and minimum values of 
E everywhere will not be too large. Statistically speaking the 
dispersion in E will not be enormous. Similarly a realistic (non 
ideal) manufacturing proces will incur variation in the depth D and 


width B of the bar (refer to Fig. 1.1). Thus B and D are to be taken 

as stochastic processes: B(z,y) and D(z,x) , respectively. Most likely 

the departures AB and AD of the width and depth from corresponding 

mean values B Q and D 0 [B = B q + AB similarly D = D q + AD] will be 

confined within a few percentage points. Now the estimation of 

probability distribution function F^in terms of F~, F^, F^ and F^will 

e c B D E 

not be a simple algebraic task. In fact the deterministic equations 
(1.1) may not even be valid for the mean values, i.e. 


? fo 

° E o 


( 1 . 2 ) 


The stochasticity in material properties and in geometry which modify 
the system stiffness is of principal importance here in the estimatior 
of randomness for strains, displacements, and any relevant response 
quanitity. Published papers, which deal with the estimation of mean 
and standard deviations (correlation matrix in the case of correlated 
stochastic variables) are reviewed in this report. The finite element 
methodology has been the focus in recommending practical solution 
strategies which consider randomness and particularly the spatial 
variability of stochastic processes in structural systems. 
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2 . Problem Statement 


We express a finite element form for the equation of motion 
of a structure in the symbolic form: 

S R « P (2.1) 

where: S: system stiffness operator 

R: system response history 

F: forcing function 

For a vector of random processes which define the system 
S, the corresponding stochastical finite element system of 
equation will then be 


S R = F 


( 2 . 2 ) 


In a generalized probabilistic finite element problem we shall 
have: 


F : the probability distribution function for basic 

X 

uncertain quantities are prescribed. 


F : the probability distribution functions for response 

R 

quantities are to be calculated. 
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From a design-analysis view point the probability that a 
certain response quantity R^ would exceed a predetermined 
prescribed value R^*, i.e. 

F ( R . * ) ; Exceedence of R: with respect to R: * 

R i 

is also very important. 

The key issues are then: 

i) Construction of the system stochastic finite element 
matrices to decribe S in (2.2) 

ii) Solution of R from (2.2) 

iii) Evaluation of the probability distribution 

function F from the solution. 

R 

In this report the effect of randomness for responses (say 

computation of F ) due to system stochasticity will be 

R 

highlighted. Thus for the majority of the problems reviewed 
herein we shall specialize the fully probabilistic equation (2.2) 
with deterministic forces: 


S R = F (2.3) 

o 

It may be remarked that the effects of random loading with 
deterministic systems, such that: 

S Q R * f (2.4) 


36 



can be found in various publications on random vibrations. 
Research with system stochasticity refer to equation (2.3) are 
rather scarce compared to quite a rich literature available on 
random excitation. From the practical consideration in the 
structural engineering application it is adequate to postulate 
narrow band processes for a stochastic quantity 5^. In some 
sense of a norm 1 . ®ne can than write: 

I A*. I! 

T-&1 1 (2 - 5 > 

1 

[A norm 1 . I is a quantification where 8X^1 is a real positive 
(non-negative) number associated with a physical variable 3 
Consequently, it is quite pertinent to propose that the system 
components and the responses as well will obey inequalitites 
similar to equation (2.5). This naturally makes the pertubation 
method very attractive in analyzing the system equation (2.3). 
From its incipience the stochastical finite element method 
resorted to perturbation expansion in forms of Taylor series 
about mean inputs in order to yield required statistics related 
to F . 

a* 

R 

In the sections that follows noteworthy papers which deal 
with system stochasticity, refer to equation (2.3), will be 
summarized with brief description of solution procedures and 
published numerical results. 
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3 . Papers of Historical Importance 

Boyce in 1961, published a paper of column buckling with 
stochastical initial displacement. ("Buckling of a Column with 
Random Initial Displacement", Journal of Aero. Sci. ). The 
eigenvalue problem related to free vibration of structures with 
stochastical mass and stiffness matrices was completed in 1968 by 
Collins and Thomson ( "The Eigenvalue Problem for Structural 
Systems with Stochastical Properties", AIAA Journal ). The 
pertubation method introduced by Collins and Thomson was later 
adopted by many researchers, such as Nakagiri and Hisada [N-4.1 
to N-4.8]. There latter authors derived the mass and stiffness 
matrices by employing the finite element method. The 
aforementioned two papers and [C-3.13 [B-3.13 are of historical 
significance in the research of structural mechanics problems 
with system stochasticity. 

From the standpoint of Probabilistic Structural Analysis 
Method (PSAM) the first paper the reviewer found of direct 
interest is by Mak, and Kelsey, in 1971 titled: "Statistical 
Aspects in the Analysis of Structures with Random 
Imperfunctions " . Cambou in 1972, employed a direct finite 
element formulation in first order stochastic analysis for linear 
elasticity problems [C-3.23* 

Mak and Kelsey [M-3.13 considered the out-of-plane buckling 
bifurcation of a column due to uncertainty in the initial stress 
distribution. This is the first published mathematical 
development with numerical results for any structural problem 
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with probabilistic consideration for the system stiffness 
matrix. The authors solved the eigenvalue problem asociated with 
the buckling problem: 

CK + E K e + II V 0 " ° <3 - 1) 

where K: elastic stiffness matrix 

K e : initial stress effect 

Kg: geometrical stiffness matrix 

and S and P are the stochastical initial force and the resulting 
buckling load. A similar development was adapted by Nakagiri and 
Hisada in their paper [N-4.4]. The detail of algebraic steps are 
furnished in section 4. Mak and kelsey in [M-3.1] presented a 
graph showing the probability distribution of failure by buckling 
and the effect of the standard deviation of the lack of fit for 
members on the probability of failure. The treatment in the 
paper are very clear and structural designer will find it 
suitable for application in practical problems. 

The Monte Carlo simulation technique with a finite element 
formulation was employed by Astill, Nosseir and Sinozuka as early 
as in 1971, refer to [A-9.1] and section 9 for detail. The 
authors devised a "front end" statistical package to generate a 
population of constitutive properties. The problem of failure of 
a concrete cylinder was considered under impact loading. Very 
encouraging results from that transient dynamic problem with one 
hundred realizations was reported. 
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4 . 


General Procedures for Stochastic Finite Elements 


In the existing literature, stochastical analysis of 
engineering systems are confined to the first- and second-order, 
second-moment approximations. The method calls for the first- 
and second-order Taylor expansion of any generic response 
quantity in terms of system random variables around the mean 
argument. Subsequently the mean and standard deviation of the 
response function in question can be estimated. This procedure 
is known as the delta method by statisticians. Hitherto emphasis 
has been placed on reliability analysis whereby the exceedance 
coefficients are estimated on the basis of means and dispersions 
of response quantities. A more precise estimation of exceedance 
calculation will necessitate the knowledge of higher-order 
moments. The existing literature on stochastic finite elements 
is rather deficient in evaluating these higher-order moments. 

The methodology for the first-order second-moment 
approximation is quite complete for linear systems. The second- 
order perturbation formulations are rather recent. Even though 
the methodology is straightforward and conceptually amenable to 
nonlinear dynamic systems, the details of computational strate- 
gies suitable for finite element systems especially with a large 
number of stochastic parameters cannot be found in existing 
literature. 

A thorough review of published research in the are of 
probabilistic analysis for finite element systems reveals two 
major directions. Theoretically, the perturbation formulation 
and numerically the Monte Carlo simulation are the only courses 
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available so far. In this section stochastic finite element 
formulations according to the perturbation method is detailed. 

The notion of finite element spatial discretization in a 
stochastic model on the basis of the scale of fluctuation is also 
reviewed here. The Monte Carlo simulation technique is more of a 
statistical method hence it is described in the next chapter. 


Perturbation Method 

The systematic development of the stochastical finite 
element formulation according to the perturbation method was 
initiated by Nakagiri and Hisada, refer to [N-4.1] - [N-4.8], 

They essentially employed the perturbation method [B-1.4] and 
reatined up to second-order terms. In order to focus on the 
stochasticity of the system, the load vector (the right-hand side 
of the equation of equilibrium) was taken to be deterministic. 

In this review, the equations furnished by Nakagiri and Hisada 
will be rewritten using the notations that appear in [B-4.1]. in 
the interest of clarity, indicial notation will be employed 
whenever required. 

The general discussion may be started by examining the 
stochastic static (global) stiffeners matrix K as an offset 
by AK from a preselected deterministic value Kq, then 

K * K Q + AK (4.!) 

Now for each element ij, the equation reduces to: 




+ AK. . 
IJ 
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A superscript will be used here to indicate the corresponding 
variable pertaining to an element "s"; hence 


cr( s ) _ ~( s ) 
K ij - K 0i j 




(4.2) 


The random variables which govern the system stochasticity , are 
collected as a vector {X} with components X^. Conceptually, both 
the global and element stiffness matrices, K and K^ s) , 
respectively, can be Taylor expanded about a preselected vector 
Xq where 


x i * X 01 + 1X i 


(4.3) 


leading to 


or 


aK i 


K ij * K 0ij + iK ij ■ X 0ij * ^ 


AX 


1 5Jx ij 

+ 77 r—b ix i ix m 

ax,ax„ 

x m 

X ij = K -ij + °ijJL AX Jt + P ijtm AX X AX m 


(4.4) 


(4.5) 


where 


5Ki 

a = 

°> l ax. 


and S,. 


"ill 


ax ax, 

in x 


It should be noted that in the above equation, the 
expressions bevond the quadratic terms are truncated by Nakagiri 
and Hisada. Ther is no such restriction (refer to Eq. 3 in 
[B-1.4] ) in a general perturbation technique. A similar 
expansion, consistent with the second-order perturbation of the 
stiffness matrix, can be implemented for the displacement vector 
U: 


au. 


a 2 u. 


0- * u rt . + — - ax. + — ax. ax, (4.6) 

1 01 J jkv 3 * 


ax. 


ax.a Xjl 
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Now a static finite element system, with a deterministic load Fq 
can be solved when the stiffness matrix and consequently the 
displacement vector are stochastical in nature. The governing 
equation of equilibrium then becomes 


K ij U j * F 0i 


(4.7) 


Now, substitution of Eqs. 4.6 and 4.5 into the above equation 
leads to 


t K 0ij * + • 




5 2 U. 


Kj + T^r 1 “x + * F 0i 


(4.8) 


&X i & V X m 


One compares the zero-th, first and second degree terms 

containing AX., AX_ , 
x m 

etc. and obtains the following recursive set of equations: 


K 0ij U 0j - F 0i 


au . 

K ° ij 8X t ' ' “ii 1 " 0 ! 

d 2 U. 

K 0ii - - 1 - t p ijto u 

J ax.ax J 

x m 


_ . + a . . 

0 j 13m 


dU. 

-rr 1 ] 

ax. 


(4.9a) 


(4.9b) 


(4.9c) 


It is interesting to note that the above system of equations 
can be solved once the Kq matrix is "inverted." Nakagiri and 
Hisada remarked that numerical calculation will be faster in 
their method as compared to a Monte Carlo simulation since the 
latter necessitates a separate inversion at each numerical 
realization for the random vector X. 
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au. a 2 u. 

Once Upi' — ”■ and are obtained by solving a set of 

dX. 3X.ax 

linear systems of equations TEq. 4.9), one can compute the 
expected value of U. The expected value operator E, when applied 
on the Taylor expanded form for U (Eq. 4.6), one obtains 


au. a 2 u. 

:|U ± ] = U oi + E [LX.) + Z E[AX. AX*] 

* V J *Y JKY J 


(4.10) 


ax . 

j 


ax . ax. 

3 i 


The authors suggested that in a "deterministic" computation with 
Kg , the stochastic vector Xg should be chosen to be the mean 
of X. Then 


E [ x i ] = X 0 . i.e. E[AX.] = 0 


(4.11) 


This would simplify the expression for the mean U in Eq. 2.10 
leading to 


a 2 u. 

E[Ui) - U 0i + - l E[ AX* AX,] 

AX AX, J 


(4.12) 


It is convenient to introduce the covariant matrix Cov[X,X] such 
that 


E[AX. AX^] = Cov [X,X] i j * Cov [X i ,X ^ ) 


(4.13) 


Now the mean displacement can be computed from: 

a 2 u. 

Etu.) * U 0i + — — Cov [X • ,X . ] (4.14) 

1 ax. ax J 

J ^ 

The second-order Taylor expansion of K and U in Eos. 4.4 and 
4.6 limits up to second-moment terms in the above equation. 
Consistent with these second-moment terms in Eq. 4.14, one 
evaluates the dispersion of in terms of the variance operator 
Var [ U i ] : 


Var [U^ = E [ U 2 ] - { E C LJ i ] } 


(4.15) 
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The second term of the right-hand side in the above equation was 
obtained from Eq. 4.14. The first term on the right-hand side is 
calcualted retaining only the second-order terms leading to 


BtU i J * U 0i + 2U 0i 


a 2 u. 

Cov[X X ] 

ax ax, 3 


a 2 u i a 2 u. 

+ 7^ cov[x. ,x,] (4.16) 

5x j 3 * 

Once the expected value E[U^j and dispersion var[U^] are obtained 
from Eqs. 2.14 and 2.16, the corresponding statistics of the 
strain e and stress t can be obtained by utilizing the strain- 
displacement transformation B in terms of the shape 
functions N and the constitutive tensor C. The algebra is sum- 
marized in Eqs. 15-21 in [B-1.4J. 

The aforementioned general technique, described by equations 
4.1 through 4.16, was illustrated by Nakagiri and Hisada in their 
first paper [N-4.1] where only the variation of the shape func- 
9 ns were considered. For a triangular meshing, a shape 
function N was written in terms of the area coordinates , L 2 
and L 3 and the nodal point coordinates. This is a standard 
finite element procedure and the details can be obtained from 
[Z-l.l], in this first paper, the stochast ici ty of the nodal 
coordinates were considered. An element stiffness matrix K (s) in 
an isoparametric formulation was obtained from the corresponding 
stochastic strain-displacement transformation matrix B^ s *, the 
constitutive matrix (stress-strain relationship) C (s) as well as 
the Jacobian transformation j (s) whose stochasticity is due to 
those of the nodal points. Integration over the element in terms 
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of the area coordinates and L 2 yielded: 

K<S> * // ig(*)]T c (s) B (S) |3r <s) IdLjdL. 


(4.17) 


where the determinant of is indicated by |j* s *|. For a 

nodal point coordinate (x,y), the Jacobian assumed the form: 


J = 


dX_ 

dL, 


ax 

aL. 


ax 

aL. 


ax 

aL. 


'2 w ~3 

which was written as 


(£L- - &L.) 

' aLj aL 3 ; 


flSL - 
^aL 2 aL 3 J 


(4.18) 


J = J Q + AJ 


s ) 

The terms in the B matrix involved expressions 
like and which were obtained as 


r 

dN 

ax 


dN 

ay 


[ J (,) ]" 1 


. ^ 
dN 

yw *S 

dN 

aL x 

dL 3 

A/ 


dN 

dN 

bL x 

- dL^ 


(4.19) 


(4.20) 


It was then possible to evaluate the and P 0 j An , terms in Eq. 

dN 


4.5, once the second-order Taylor expansion of |J|, were 
obtained. The authors denoted 


where 


+ + Dj (4.21) 


D 1 * J ll iJ 22 “ J 12 J 12 AJ 21 ” J 21 AJ 12 + J 22 AJ 11 (4.22a) 

D 2 * AJ n AJ 22 ” A J 1 2 A J 2 1 (2.22b) 
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e ) 

The required inversion of J v in Eq. 4.10 was expressed as 


* (1 ' Wo 


IJ]* 1 = (1 - fj^|- - p^T + Tj^ T t ' 


A j 22 " AJ 12 


" AJ 21 AJ 11 


J 0 22 ** J 0 1 2 


_J 021 J 011 


(4.23) 


Finally a... g . . . tensors were obtained after Taylor expansion 
1 ljt ljAm 

of x and y up to second-order terms. 

An example was illustrated where the nodal coordinates were 
taken as stochastic processes defined by a power spectrum. A 
homogeneous Wiener-Khintchine relation was assumed for the 
correlations Cov[x,x], Cov[x,y], Cov[y,y]. 

The autocorrelation r(|x£ - Xj|)for a homogeneous stochasticity 
was obtained in the following form: 


R(|x^ - x j | ) =2 f s(X)cos 2X.|x^ - x j | dx (4.24) 

from a given spectrum s(X). 

The paper does not present detailed numerical results. The 

computational procedure for the a..„ and 6 . . . tensors are not 

i j x i j xm 

discussed either. It should be noted that in a practical finite 

element formulation with stochastic variables the computation of 

a. and g. . would demand substantial numerical and programming 
ljJt l jem 

effort. 
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In their second paper [N-4.2] Nakagiri and Hisada considered 
stochasticity in the static stiffness matrix K due to 

(i) variation of constitutive properties, where C [ X J is 
considered 

and 

(ii) variation in boundary data. 

The general methodology described before is implemented for those 
two cases. The paper details out plane stress/strain examples. 
These steps are crucial in developing a stochastic finite element 
code with plane elements. However, proper adaptation of the 
algebraic derivations to general finite element stiffness 
matrices (and to mass matrices as well) could lead to the 
formulation pertaining to three-dimensional solid and plate or 
shell elements. In the interest of focusing on the method the 
two-dimensional linear elasticity example will be sketched out 
here . 

The stochasticity of the constitutive properties was 
considered first. For a plane stress/strain element the Young's 
modulus E and the Poisson’s ratio v were introduced as bivariate 
stochastical processes, in the form of E(X) and v(X). The random 
vector X is indeed dependent upon the spatial coordinates x^ and 

x 2 * 

The element stiffness matrix is composed of 2x2 submatrices 
obtained from two shape functions (x^) an< ^ re ^ er to 

[Z-l.l] for details. This submatrix can be written in the 
following form: 
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aN . 

dN . 

aN i aN . 

i 

— — 1 + V ' 


aXj 

ax : 

axfj ax 2 



(4.25) 


aN. 

8N . 

aN i aN^ 


.axj + 

axy 5y 2 


Symmetric 

5N^ aN,. aN^^ SN ,. 
v ' 5x 1 + TyT 2 dXj 


The values of the stochatic variables X, v' and v" are expressed 
in terms of E and v for the plane stress/strain cases: 


variable 


plane stress 


plane strain 





E ( 1-v) 

(1 + v) ( 1-2 v ) 


(4.26) 

1 - 2 v 

2 (1-v) 


V V v 

1-v 

As before the Young's modulus and the Poisson's ratio is Taylor 
expanded about their means only up to the first order terms: 


E 



(=E q + AE) 


v n + AX. 

0 3X. 1 


(4.27a) 


(4.27b) 


This leads to the following form for the element stiffness 
matrix, K^, with second order terms: 
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s(s> , K (s) + SlU t AE (S> 
0 


* Mil’ iv (Sl 


»~( s ) 

a v 


*2 K (S) 

+ 5E (Sl ^ lS ’ 


(s) »„<•> + 1 iLilfi (S v (s ) 


AE Av 


2 & v (5)2 


(4.27c) 


[No sum over repeated index "s"] 


It is to be noted that the stiffness matrix K is proportional to 

5 2 K . 

the Young's modulus E hence jv— is zero. 

5E 

The displacement vector U when perturbed up to second order m 
terms of AE and Av leads to 


U = U n + 


au 


0 dE (S) 


AE 


(s) + 5U 


av (s) 


Av 


(S) 


9 

k *• 


(S) AP (t) 


+ l [ , 5 , U 7 7 T AE'*" * E ' 

2 5E (S) 6E (t) 


a 2 u .„(s) A Jt) 


(4.28) 


+ 5E (S) av (t) 


AE' ' Av 


5 Av (s) Av (t) ] 


+ av (s) av (t) 


The expansion for the displacement vector involves summation over 
all elements (as described by the superscript "s" and "t ), 
whereas that for K (s) pertaining to a particular element s is 
described by variations of E and v in that region. 

In the case of a deterministic load vector Fq the unknown 
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partial derivatives of U with respect to E and v can be obtained 
by considering terms with AE and Av in 


K U * F 


0 


leading to: 


U 


0 



F 


0 


5U 

BE 



^O 1 ' 1 

[ Kq ]- 1 

IV ' 1 



rlE Mi 

BE BE 


B 2 U _ r „ ,-1 , BK BU . BK BU 

~ ~ " “ l *n J 1 ~ ~ “ ~~Z 

BE 3 v BE Bv Bv BE 


+ U Q ] 

BEBv 





U, 


BK BU j 
Bv a v 


(4.29) 

(4.30a) 
(4.30b) 
(4.30c) 
( 4 . 30d ) 

(4. 30e) 
( 4 . 3 0 f ) 


The authors suggest the computation of the mean and dispersion 
of U from the above expressions. As claimed by the authors to be 
a strong point of their formulation the aforementioned equations 
involve the inversion of a single [Kq] matrix. (In a Monte Carlo 
simulation each realization would demand a separate inversion 
of K. ) 
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The strain-displacement transformation matrix B is 

. ( s ) . 

deterministic hence the strain vector e becomes: 

c (s) = b ( s) U (4.31a) 

= B (s) ^ . . . as in (4.28) ] (4.31b) 

0 5E 

hence the mean strain E [e] and its dispersion Var [e] can be 
calculated directly. Finally the stress calculation involves the 
stochastic constitutive matrix C: 


C 


= l 


1 

v 


0 


1 

0 



(4.32) 


Employing the explicit definition of C from (4.26) one obtains 

•k p * p 

directly those partial derivatives like — and — . Substitution 

5E &v 

of these Quantities lead to the expression of mean and dispersion 
of the stress components. 

The authors have not commented on the numerical 
implementation of mean and dispersion calculation of the stress 
vector a = C B U (4.33) 

The rest of the paper [N-4.2] elaborates the concept of 
adopting a stochastical description of the boundary data. The 
nodal degrees-of-f reedom (with the prescribed stochast icity ) 
which pertain to the boundary were designated with a 
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superscript ^ and the remaining (interior) degrees-of freedom by 
1. Then the static stiffness matrix, the displacement and the 
load vectors could be partitioned leading to 


\ (11) 

K< 12 ' 


r«»>i 



< 

t 

.it' 21 ’ 

k (22 » 


[u' 2 >J 



(4.34) 


Thus the unknown displacement vector (associated with the 

~( i ) 

interior nodes) U becomes: 


~( 1 ) 
U 


[ K ( 1 1 > ]“ 1 [F (1) - ic< 12) U (2> ] 


(4.35) 


The authors pointed out that the aforementioned equation 

^ / 1 \ 2 ) 

indicated a linear relation between U and U hence 

conjectured the possibility of numerical computation of the mean 

~( 1 ) 

and dispersion of U from those of the right hand side 
quantities from equation (4.35). Finally relevant statistics for 
the strain e and stress o distributions could be obtained 
according to the equations (4.31) through (4.33). 

The authors do not include specific numerical examples for 
this problem of stochastici ty with random boundary data. As in 
[N-4.1] and in (4.24) a power spectral density function in the 
Wiener-Khintchine form was suggested to account for the spatial 
variability of the stochastic quantities. The authors did not 
elaborate on numerical computations of means and dispersions of 
the stress components from those of the given constitutive matrix 
and calculated displacement vector, refer to (4.33). 
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The possibility of extending the perturbation technique 
sketched out in [N-4.1] and [N-4.2] for nonlinear problems is 
discussed by Nakagiri and Hisada in their third note. A specific 
nonlinear constitutive model in the following form was explored: 


o = fU> = E e -E b[n{ 1 + (^)} a -l] 

b 


(4.36) 


A nonlinear constitutive tensor C could then be assumed in terms 
of the effective strain e and effective stress o, which are 
defined to be 



(4.37a) 


and 


o = f (e) 


(4.37b) 


~ * 

In principle , for a selected value of U to be D the 
constitutive matrix C (s) for an element "s" was Taylor expanded 
up to second order terms with respect to the stochastical 
constitutive variables a and b in (4.36) in the following form: 



(U ) 



( U* ) + 



a a 


* Aa 

U 
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ac 

5b 


s ) 


* Ab 


u 


, .2 ~(s) . -2 s(s) _ _ 

+ I [- — ^ — Aa 2 + 2 - — + Aa Ab 


2 

a u 


aa ab 


* 2 s(s) , 

* A-r§ — (4b) 2 ] 

ab^ 


(4.38) 


^ ( s ) * 

The stochastic element stiffness matrix K was defined as a 
quadrature of [B* S M T [C (s) ] ( B ^ s ^ ] at selected Gauss integration 

point leading to a form: 


K<s> , K (S» + » g'*j *,(•) + Mill ib < = ) 

0 ~ ( S ) - r-( S ) 


a a 


♦ I [^1 Ua <*>> 2 ♦ 2 

2 aa (s) 


ab 


(4.39) 


a 2 k (5) 

Sa ls> 5b (s) 


4a (s) 4b <s) 


jl 2 ar( S ) . . 0 

o K r Awls) ^2*1 

’ ] 


In their derivatives of the stiffness matrix could be computed in . 
the following form: 


8K (S) 

aS (s) 


- / ib (s) ] t [ --£ - !v! ] i® 1 dv 


a a 


(4.40) 


*2 ~( s ) . . _ 2~(s) 

i rnr ■ 1 IB 1 [ 4 si IB1 dv 

a a 1 ; a b ts ' v a a lsJ a b (S) 


Thus all partial derivative of the global stiffness matrix K can 
be obtained by assembling the aforementioned corresponding 
partial derivatives defined for each element. Along with the 
Taylor expanded version for K-- in the same form (4.39) — second 
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order perturbation of U in terms of U 0 , 22 and 22 leads to the 

da 6b 

following system of linear eauations: 


K° 

* 0 

(U °) 

U J = 

F 

K° 

* 0 

(U u ) 

22 + 

5K 



5a 

5a 

K° 

(U*°) 

22 + 

*** 

5K 



6b 

6b 

,2 

_I u° 

+ &£ 

5U 

6a 

6b 

da 

5b 


4. 

K° -2! 

*** 

U 



5a 

5b 


(4.41a) 


(4.41b) 


(4.41c ) 


( 4 . 4 Id ) 


This permits computation of mean and dispersion of U in terms of 

0 *%■> 

& U & U 

~~Z' ' ~Z T. an< ^ i- n terms of means and dispersions of a and E. 

6b 5a 5b 

Numerical evaluation of the derivatives, according to (4.40) 
could be somewhat complicated when the nonlinear stress-strain 
relation (4.37b) is elastoplastic in character as in (4.36). 

1 his step will consume substantial computational resources in a 
large problem. 

The paper does not elaborate on numerical implementation of 
the steps presented therein. 
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Nakagiri and Hisada applied the perturbation method to 
evaluate safety and reliability for finite element representation 
of structural systems. The paper [N-4.5] described a framework 
to apply the mean and dispersion of response quantities, as 
evaluated in [N-4.1] through [N-4.4] , in order to calculate 
safety indices. The methodology of the standard reliability 
technique applied therein can be found in [W-4.1]. 

The time history analysis with a stochastical description of 
a proportional damping matrix was presented in [N-6] . Some of 
the crucial aspects of the latter are described below. 

The equation of motion for a damped finite element system 
can be written to be 

MU+KU+K' U » F ( t ) (4.42) 

In the case of proportional damping the damping matrix K' is a 
linear combination of the mass and stiffness matrices M, K, in 
the form: 

K' = a M + b K (4.43) 

In the specific example [N-4.7] the authors considered determi- 
nistic mass and stiffness matrices, then M = M and K = K and 
focused attention on stochasticity of the damping matrix via the 
random variables a and b in (4.43). In the computational step 
that was presented in [N-4.7] the authors formulated a broader 
class of problems where the damping matrix K* was decomposed into 


★ 

conventionally C is used to indicate the damping matrix. In 
order to avoid confusion with using C for the constitutive matrix 
in t B— 1.4], a nonstandard notation, K' is used to denote the 
damping matrix herein. 
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a deterministic component K ' q and a stochastic part which is 
proportional damping in nature i.e.: 

K' = K' 0 + aM + bK 

Furthermore, K'q was not necessarily in the form of a Cangley 
series with K and M. Hence Kg was not reducible to a diagonal 
matrix with real mode shapes pertaining to those of the 

undamped system. A generalized version of the aforementioned 
equation would be 

K' = K' 0 + I a i K' i (4.44) 

where each K'. would reduce to a diagonal when transformed into 
the modal coordinates as follows: 

[ 4] T (K 1 i ] U] = <k[> (4.45) 

n 

where <K^> is a diagonal matrix. It may be remarked that such a 
generalization is very appropriate for a wide class of practical 
problems. 

In the formulation a generalized coordinate {q^} was Taylor 
expanded in terms of the stochastic parameters a^ of (4.44) up to 
quadratic terms: 

„ Mq^} i 

<<V * {q 0 i } + — =*- Aa i + i 1 ~ ■ V < Aa i> (4 - 46) 

aa j Sa j da X 
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substitution of (4.46) and (4.45) in (4.42) led to the following 
set of equations when the modal representation was sought: 

m <?i + <>i kg Qi + k q. = 4>T f (t) (4.47a) 


m 


i 





i da . 
1 


k'! 

3 


(4.47b) 


m 


da . da „ 
3 1 


T 

♦i k o 


s2 qj 

'i da . da . 
3 * 


(4.47c) 


+ k 


a 2 ., 


e V a * 


dq. dq. 

- (k" + k" — i) 

1 da i &a.' 

* J 


(note: no sum over repeated indices) 

The modal mass and stiffness components were obtained from the 
mode shapes as: 

♦i M +i * M i and ♦[ K ♦i - K. (4.48) 

The authors solved the aforementioned set of equations by 
employing Newmarks implicit time integration scheme ("p * A"] . 
Numerical results for a tower with fourteen beam elements 
subjected to El Centro (1940) NS acceleration input was selected 
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to be the input ground motion. 

The statistical computation was simplified by assuming the 
random variables a in (4.44) describing the random coefficient 
for the proportional damping matrix to be of zero mean Gaussian 
distribution. Thus the required input were Cov [a^,aj]. It 
should be noted that the third moment: 

E [a ^ ,a ^ , a ^] = 0 (4.49) 

and the fourth moment reduced as: 

E [a ^ ,a j ,a^ = Cov [a^a^] Cov [a Jt ,a m ] + (4.50) 

Cov [a. ,a^) Cov [a^,a m ] + Cov Cov [a.,a A ] 

Numerical results in the form of graphs indicated expectation and 
(» 3 _ a ") bounds of top deflection and the effects of Cov [a^,aj] 
on the standard deviation of top deflection for the tower 
problem. These are perhaps the only meaningful numerical results 
published for dynamic analysis of a finite element system with 
stochastic damping matrix. 

A column buckling problem [N-4.3] with stochastic 
description of the stiffness matrix K and the geometrical 
stiffness matrix led to the computation of the buckling load 
via the following eigenvalue problem: 

[it) (?) - X [Kg] (?) (4.51) 

The buckling load is related to the eigenvalue \ and the bent 
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shape is described by the eigenfunction {z } . It should be 
remarked that the content of this paper [N-4.4] is identical with 
[C-3.1] where the free vibration problem: 

IK) { z } = w [M] { 2 } 

was described. In (4.52) w and {z} are the stochastical natural 
frequency and the mode shape due to stochastic mass and stiffness 
matrices, M and K, respectively. 

In the note [N-4.4] the authors presented the problem of 
buckling of a cantilever beam with the stochastical descriptions 
of end restraints as shown in the Fig. 4.1: 



Fig. 4.1 Cantilever Beam Buckling Problem 


The element stiffness matrix for element number 1 which was 

attached to a stochastical spring with translational spring 

s El **' 

constant ) -r and rotational spring constant (— £— ) would 

1-s L J l~c L 


be 
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In the procedure that followed the element stiffness matrices 
were assembled in global matrices [K] and [K g ] . The stochastic 
processes, viz [K] , X, and (z) were Taylor expanded up to second 
order terms with respect to the random parameters s and c about 


their mean 

values. 

Sq and Cq, respectively. Thus 


0 

K « K 


+ ™ 4S + i l + U =) 2 2 

as ac asac 

( AS) ( AC) 

+ Ac 

be 



+ ^s) 2 ] 

a s* 

(4.55) 
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since K in (4.53), had rather simple algebraic expressions in 

terms of s and c evaluation of the partial derivatives, 

i.e. a, 3 tensors of [B-1.4] , are indeed straightforward. Thus 


dK _ El S 0 

~ = “3 

as L (1-S Q 


2~ S 

S K EI s o 

5 T 

Ss L° (1 -Sq) 


5K _ EI C 0 

~ = ~ 2 

ac L (1-C 0 )^ 


2— C 

a K _ EI ^0 

~2 ~ ~ T 3 

ac (1-C 0 r 


1 0 
0 

Symmetric 


0 

0 

0 


1 0 
0 

Symmetric 


0 

0 

0 


0 


0 


0 

0 


Symmetric 0 


0 


0 


0 


L 
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(4.56a) 


0 

0 

0 

0_ 

(4.56b) 


0 

0 

0_ 

(4.56c) 


0 

0 

0 

0_ 

( 4 . 56d ) 


Now similar expansion, as in (4.55), for the eigenvalue X (which 

~ 2 

was proportional to the buckling load P, P = and the bent 
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shape { 2 } were carried out in the following form: 


x n + is + £ 4C * 4 4| “ 2 

0 65 sS 2 SC 2 


+ ill is 4= ♦ i 4C 2 

iHa? 2 &C 2 


(4.57) 


( is ) 2 + 


— (AS) (AC) ^ (^) 2 

»sa? 2 ac 2 


(4.58) 


2 ~ 2 ~ 2 

Note eventhough - 0, and are nonzero due to 

asdc asac asdc 

coupling through the implicit "inversion" in a linear eigenvalue 
problem. Substitution of (4.53) through (4.58) in (4.51) led to 
the following system of linear equations (when the coefficients 
of As, Ac, (As) 2 , (Ac) 2 are set to zero individually) 


( [K 0 ] - X 0 [K g ]) {Z 0 } = 0 


(4.59a) 


aK 

ax 

*** 

as 

as 9 

5K 

ax K 

ac 

a? 9 

2 ** 

2~ 

a k 

ax 

~2 


as** 

as* 


SK ax 

<S» 'W 

as as 


(4.59b) 


(4.59c) 


( 4 . 59d ) 


+ (K 0 - ^0 K g> ( ^! } ’ 0 
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( 4 . 59e ) 



( 4 . 59f ) 


It is to be noted that computation of partial derivatives such as 

~~z~~Z et c. in the aforementioned equation, can be carried out 
ds 5s3c 

by solving the generalized version of the linear eigenvalue 
problems of the form: 


{U} + [A] {z} + {V} = 0 (4.60) 

where X (scalar) and {z} (vector) are unknowns but { u } and {V} 
(vectors) and [A] (matrix) are prescribed. 

The authors presented numerical results for a sample case 
and demonstrated the accuracy of this second order perturbation 
method. 

Finally, the mean the dispersion of X and {z} were obtained 
following the methods in the previous paper [N-4.1], (N-4.2J and 
[N-4.3] . 

Nakagiri and Hisada also employed the perturbation technique 
to the specific cases of stochastic Winkler foundation [N— 4.6] 
and for random misfit in frame structures [N-4.8J. These two 
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papers are of peripheral importance in describing the general 
procedure of the stochastic finite element formulation with 
perturbation techniques. 

In a conference paper (N-4.9] the authors summarized the 
perturbation technique developed in [N-4.1J to (N-4.8]. 

Numerical results for two specific problems were presented. The 
expected value and dispersion of stress intensity factor for an 
edge crack with uncertain length were presented in graphical 
form. The second problem dealt with mean and standard deviation 
of inplane stress developed in a long strip. The authors 
compared the results of the first and second order approximations 
(where the stochastic processes were Taylor expanded up to linear 
and quadratic terms, respectively). In certain cases the 
difference of result was quite significant. Hence the authors 
recommended the formulation with second order approximations. 
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Handa [H-4.1] initiated stochastical finite element 
calculations to enhance design analysis of civil engineering 
structures. In a sequence of research reports and conference 
papers [H-4.2], [H-4.3], [H-4.4] Handa and his associates 
employed finite element analysis technique to estimate expected 
values and correlation coefficients of static stresses and 
displacements for trusses, frames and beam structures. 
Stochastical variations of structural section geometry, material 
property as well as that of applied loading were considered. All 
their discussions were restricted lognormal distribution of 
stochastical parameters. For example, for a finite element (s) 
the carrying capacity R (s) and any load effect S (s) were assumed 
to be lognormally distributed. This facilitated the construction 
of the safety margin z defined to be: 

z (s) = In R (s) - In S (s) (4.61) 

to have a normal distribution. 

The presentation [H-4.1] detailed out the first order 
perturbation method. The authors remarked that the error 
associated with neglecting the higher order terms will not exceed 
20% at most. In the interest of brevity the steps are not 
repeated here since more detail algebraic development are 
presented in this section in equations (4.1) through (4.16). 

The authors presented several numerical examples of the 
first order second moment formulation. Two noteworthy cases 
among those will be summarized here. 
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Fig. 4. 2 Cantilever Beam Problem 

A cantilever steel beam, as shown in Fig. 4.2, with deterministic 
length L and diameter D, was analyzed with the following 
uncorrelated stochastical processes: 

( i) loading W(x) , 

(ii) second moment of area I(x) and 

(iii) modulus of elasticity E(x). 

The mean values Wq, Iq, Eq and the standard deviations 
V a i and a E were P rescr> ibed as input data. The spatial 
variability of the aforementioned stochastical processes was 
taken to be exponential. The autocorrelation functions for such 
processes were expressed in the form: 

P (x^Xj) = exp (-k | x i - Xj|) (4.62) 

The constant k (with dimension of length inverse) was anticipated 
be obtained experimentally. Numerical computations were 
carried out by assigning k * 0 (fully correlated), K ♦ - 
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(uncorrelated) and intermediate k = 2 (partially correlated) 
cases. It was demonstrated that the standard deviation for the 
displacement of the tube at the free end was significantly 
dependent upon the spatial variability criteria as depicted by 
the correlation coefficient in Fig. 4.3. 

v* 


laticn 


Partied Correlation K = 2 
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second moment of areas and the moduli of elasticity were 
considered as random variables (with prescribed means and 
standard deviations). The stochasticity of fully correlated 
loading was expressed by a fully populated 14 x 14 matrix. 
Correlation between cross-sectional area, second moment of 
inertia and moduli of elasticity i.e. Cov[A,E], Cov[A,I], 

Cov[I,E] were taken to be zero. The variance matrix for each 
random quantity (associated with those 18 nodes) was thus 
obtained in the form of a diagonally dominated banded matrix. 

The authors used a computer program [H-4.3] to carry out the 
stochastic finite element calculations. Mean values and 
covariance coefficients of displacements of all nodes and 
stresses in each element were evaluated. The correlation matrix 
for stresses was fully populated but decayed with geometrical 
distance between two elements. 

Scale of Fluctuation Method 

Vanmarcke conducted extensive research (V-4.1], [V-4.2J, 
IV-4.3], [V-4.4], on engineering problems dealing with partial 
differential equation of mathematical physics where the 
properties of the domain are random processes in the spatial 
coordinates. The outstanding contribution was to systemically 
develop expressions for local spatial averages of stochastical 
quantities as well as variances. For example, in one— dimensional 
situation, a stochastic process is a "random function" of the 
spatial coordinate x: X^(x). In general the expected value for 
the product X^Xj) and X i (x 2 ) will be 
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(4.63) 


E [X.fXj) X.(x 2 )J 
or 


Cov [x., x j / x 2 1 

Cov [X., | x 2 ~ x j | t Xj J 


For .l^ tionar y processes (with respect to spatial coordinates) 

the values in (4.63) are independent of x x and could be expressed 
as 


< X i 


(X. ) X ; 


(x 2 )J = 


Cov 


[X., 


► ” x l I 


(4.64) 


or 


E [X. (x x ), x. (x 2 )] 


i P X t (x 2 “ x l ) 


(4.65) 


where p is the autocorrelation function. 

Within the framework of finite element analysis the 
distributions over an element "i" are "smeared out" and 
statistical average quantities were defined over the element 


domain Lj in the form 


X 


ij 


jth element 




(x) dv 


(4.66) 


Now x.. is the equivale nt random variable over jth element 
associated with the stochastical process X. (x). if x. i s a wide 
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sense stationary 
and variance Var 


process 


,x ii' 


then one defines the mean m 
in the form: 



- 


m~ 
x . 


Var 


lx u> 



(L.) 


(4.67) 


in which Y ~ (L.) described the dependence of the element size on 
x i 

variance function. From the identities of random signal process- 
ing one could write: 


x x 

Y~ (X ) = / / p~ (x -x_) dx. dx. (4.68) 

x i 0 0 i 1 z 


Furthermore, the variance function y~ (x) has the property: 

x i 


ytr (0) = 1 and at y~ (x) = 0 (4.69) 


The principal contribution of Vanmarcke's presentation is to 

define the scale of fluctuation 0~ based on the asymptotic 

x i 

behavior of y in the following form: 


>y (x) * as x >> 0~ (4.70) 

x i x x i 


and 


Q~ = 2 / p~ (x) dx 
x i 0 x i 


(4.71) 
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whenever the above limit (4.70) and the integral (4.71) exist 

An important interpretation of ©~ is in terms of the Wiener- 

x i 

Kinchine spectral density function g~ (u) where 

X i 


9y ( w ) ~ ~ j Prr (x) cos ux dx 
X i * 0 X i 


(4.72) 


and then 


9 ~ 


x 


i 


= itg~ (0) 
i 


(4.73) 


It is important to note that during the selection of mesh size, 
for a stochastic process indicated by Wiener-Kinchine spectral 
density the characteristic length of a finite element region 
should be less than the scale of fluctuation. 

Considerable development regarding the scale of fluctuation 
for general two-, three- and n— dimensional random processes 
could be found in the textbook [V-l.l], Vanmarcke also presented 
very useful approximate formulae for the scale of fluctuations 
based on the asymptotic behavior of the correlation function. 
Detail mathematical development for unidirectional and two- 
dimensional random variates along with useful algebraic 
identities could be found in the conference paper [V— 4.1], 

The journal paper [V-4.2] described in detail the problem of 
static deformation of an idealized shear beam shown in Fig. 4.5. 
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Fig. 4. 5 Shear Beam Problem 

The continuum equation for the deflection v(x) was expressed as 

dv ffi * <*(x) V(x) (4.74) 

in which V(x) was the shear force and a(x) was the shear 
flexibility. (For a rectangular beam of area A(x) and shear 
modulus G(x) the shear rigidity is given by G(x)A(x).) 

Then 

a(x) » 1/ G(x) A(x) ) (4.75) 

For the fixed end beam since v(x*0) = 0 

X X 

v ( x ) = / ~ dx * / a(x) T( x 1 ) dx* (4.76) 

0 a 0 
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Vanmarcke and Grigoriu illustrated the problem with deterministic 
T ( x )• The shear flexibility factor a (x) was assumed to be a 
stationary random process in the spatial variable x with a mean 
m<x and variance o * , autocorrelation function p 2 (x) and scale of 
fluctuation 0^. The integral relation relating the random 
displacement v(x) was given by 

V(x) = / a (X* ) T (x’ ) dx’ (4.77) 

The above equation indicates a linear transformation of T(x) into 
V(x). Thus 


E ( v( x ) ] * E[f a (x‘) T(x') dx'] 
0 


* m f T(x') dx' 
0 


and 


Var [v(x) ] = a* 


x x 


/ / P_ (X. - x' ) T (x' ) 

0 0 1 z 1 


T (x 2 ) dx| dx 2 


Some special cases of interest would be 


(4.78a) 

(4.78b) 


(4.79) 


(i) completely correlated case p~ (x) = 1 

OL 


(4.80) 


Var [ v ( x ) ] 


a Z 2 l J T (x 1 ) dx'] 2 
0 


(4.81) 
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(ii) uncorrelated case 


p~ (x) = 6 ( x ) [ 6( x ) : Dirac’s delta (4.82) 


then 


Var tv(x)) 



[T(x' )] 2 dx’ 


(4.83) 


The aforementioned evaluations of the Var [v(x)J for the general 
case (4.79) was simplified by introducting the notion of scale of 

fluctutation 9~. 

a 

Vanmarcke and Grigoriu considered the one-dimensional shear 
beam case (refer to Fig. 4.5) for the stochastic variable of 
shear flexibility a(x) subjected to two loading cases. For the 
concentrated load P, applied at the end of the cantilever beam 
the expected value and the variances were obtained in the 
following closed forms: 


E (v N ) = m~ PL 

Var (v N ) » o~ 2 P 2 L 2 (L) 

In the case of a uniformly 
Pq the authors furnished: 

x 2 

E (v(x)] * m~ p Q L (x - 


(4.84) 

(4.85) 

distributed (deterministic) load 

(4.86) 
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Var [ v ( x ) ] 


= °£ P 0 


. 2 t , 1 

L J dx. 
0 


x 


/ 

0 


(L-xjHL-Xj). 


P 


a 



x 2> d 


X 


1 

2 


(4.87) 


Numerical results in graphical forms summarize the variation 

of the standard deviation of the end displacement with L/0~, for 

a 

various autocorrelation functions such as 


p :< v x 2 > = exp I — ^©-^1 2 ) 

p-IXj.Xj) = exp (-2 

p :< x i' x 2 1 * 1 + 4 + *»p (-< l^ 2 !) 


(4.88a) 

(4.88b) 

(4.88c) 


The crucial steps in a stochastic finite element modeling 
would then be: 


Step - Is Divide a domain D into elements and 
define the element flexibility quantities as averages: 

“i 3 9 / at(x) d v / / dv (4.89) 

D i D i 

Step — 2: Calculate the mean and covariant matrices: 


{m 


a 


1 



• • • m 



m~ (1,1,1, ...) 


(4.90) 


and Cov 




£ 


(4.91) 
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For one-dimensional cases the authors used the approximate 
expression: 


o~ 

a 


a f , . , , l r ( k-1 ) Li 

Cov (a.,J.) * -y- {(k-1) J 


i 3 


(4.92) 


- 2k T, #1 - fk*!)' v a [AS2tf-±]» 


( k+1) L- 


(4.93) 


Where k * | i-j | and the beam of length L was divided into N equal 
segments . 

Step - 3: The nodal loads are to be defined by introducing 

the shape functions: 


Q. « f N. T(x) dv 
1 °i 1 


Q = {Q.} 


(4.94) 


(4.95) 


Then for the specific case of the shear beam 


v i = 2 i' 0'^} 


(4.96) 


*** 


“N } 


Then the expected end displacement and Variance of the nodal 
displacement vector: 


E [v^] ~ (2^, ••• 


(4.97) 
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(4.98) 


Cov = ° T E 0 

The approximation for the covariance matrix E will depend on o~ 
and the scale of fluctuation ©~ for a selected number of 
discretization N. 

This paper [V— 4.2] will serve as a basis to approximate the 
covariance matrix on the basis of scale of fluctuation. In 
practical finite element mesh design the scale of fluctuations 
for the random field will guide the selection of mesh size. 
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5 . Finite Element Stochasticity by Simulation 


Astill, Nosseir and Shinozuka [A-5.1] completed the problem 
of wave propagation through a random medium by employing a finite 
element modeling. Instead of resorting to the pertubation method 
to develop the stochastic mass and stiffness matrices the authors 
utilized a direct Monte Carlo simulation procedure. This could 
be the first published paper where displacement, strain and 
stress histories were generated by solving the dynamic euation of 
motion of a finite element model with stochastic parameters. The 
authors focused their attention on the impact problem and 
captured the effcts of randomness in the contitutive properties 
on the propagation of stress pulses. 

The authors modeled a concrete cylinder with 64 axisymmetric 
rings. A quadrant was modeled with 85 nodes. Uniform stress 
impact in the form of a triangular shaped pulse was considered. 
Graphs of propagation of the stress pulses were presented at 
various sections. 

Spatial variability of Young's modulus E and density p were 
considered. Numerically 100 test samples were recreated by 
employing the Monte Carlo simulation technique. Sample 
realizations of E and p were plotted against the corresponding 
mean values. Deviations from the deterministic case for the 
axial stress distribution were also displayed. The means and 
standard deviations for the octahedral shearing stress and 
miximum shearing stress were calculated using the simulated 
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population. Since these physical quantitites govern failure 
conditions in concrete cylinders the example is indeed of 
practical interest. 

The method of simulation for two one-dimensional random 
processes f±(z) and f 2 (z) (say the density and the ultimate 
strength which could vary only axially along z-axis) was 
presented in satisfactory detail. For homogeneous processes the 
cross-correlation matrix was defined as: 


E [fj^z) f^z + C)]* E Cf 1 (z) f 2 (z + C)] 
symmetric, E [f 2 (z) f 2 (z + C)] 



(C) 


°1 ° 2 r 12 (C) 


symmetric 




( 5 . 1 ) 


In there r^j(z) are the normalized auto-correlation functions. 
This matrix can be estimated from experimental data. The Wiener- 
Khinchine transform of the above correlation matrix is the 
following cross-spectral density: 


° 1 9 11 ( ^ 
Hermetian 


a l 0 2 9 22^ 

o 2 2 g 22 (T i ) 


( 5 . 2 ) 
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in which 


9 ij (»l) - 7 ; e *P dC 

Implentation of FFT (Fast Fourier Transform) based algorithm led 
to the simulated variables for homogeneous multivariate Gaussian 
processes. The density and crushing strengths were obtained by 
simulation and then the Young's modulus for each sample was 
calculated by their nonlinear transformation (as is common in 
concrete failure analysis). 

Computer code to carry out conventional finite element 
dynamic calculations was proposed whereas very sophisticated 
simulation techniques were used to generate sample finite element 
system (mass, stiffness) matrices. In order to adhere to the 
prescribed spatial distribution of stochastic processes, which 
represent randomness of material properties, the authors 
constructed cross-spectral density matrices. The required 
mathematical treatment demands thorough training in computational 
statistics. It should be remarked that merely ad hoc generation 
of realization for system matrices will prove to be completely 
useless. In structural reliability assessment the randomness of 
the system should be viewed in the light of multi-variate and 
multi-dimensional processes [S-5.1]. Gaussian processes with 
ARMA (AutoRegressive Moving Average) representation [S-5.2] are 
very useful indeed. However, for non-Gaussian stochasticity the 
computational complexity and the requirement of theoretical 
background in computational statistics could make an analysis 
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almost prohibitive at the existing level of technology. 

The aforementioned paper by Astill, Nosseir and Shinozuka 
could be the only complete treatment on simulation to be of 
practical significance. Engineers undertaking Monte Carlo 
simulation for spatially varying random processes will find that 
presentation extremely useful. It should be remarked that the 
proposed simulation technique demands advanced training in 
computational statistics especially in random process analysis. 
However, the method to generate statistics (means, dispersions, 
etc.) is straightforward once simulation techniques are 
mastered. Thus the appropriate steps will be: 

(i) to obtain a realization of geometrical and material 
properties, etc., according to design statistical 
criteria; 

(ii) to carry out conventional finite element analysis; 

(iii) to generate a population by repeating (i) and (ii) 

(iv) to construct mean, covariance matrix, skewness, etc. 
from the results of (iii). 

The mathematical treatment of Monte Carlo simulation is 
arousing new interest since the emergence of parallel processors, 
[K-7.13. Research is underway to reformulate simulated finite 
element models in order to take the advantage of inherent 
parallelism in finite element formulation [S-8.1]. 
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6 . Papers of Special Interest 

Der Kiureghian applied the finite element method to analyze 
reliability aspects of linear structures consisting of random 
variables [D-6.1 and D-6.2]. Consistant with the notion of 
computing the performance index of a structure, a stochastical 
vector s is defined to represent the effects of random load and 
random system properties. Any response quantity, like stress, 
deformation, can be included in this vector. The paper 
elaborates the first-order reliability approach, which relates 
the S vector with the allowable "strength" variables vector R 
(which typically include design stresses, tolerable deformations, 
etc.). Description of the stochastic finite element formulation, 
as applied to a linear static system, is extremely clear in the 
presentation. It may be remarked the [D-6.1] and [D-6.2] are 
perhaps the only two papers in the field of probabilistic finite 
elements, where all the conclusions and statements are 
substantiated with numeric developments. The papers are devoid 
of conjectures and ad hoc promises regarding the computability of 
large systems with probabilistic variables. The beam example 
presented in these two papers [D— 6.1] and [D— 6.2] which are 
essentially the same, is summarized below. 

A beam element in a two-dimensional configuration was 
described with three degrees-of-f reedom at each end. Associated 
with these translational and rotational deformations, the static 
stiffness matrix for the uniform section was described by the 
following stiffness matrix. 
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Fig. 6. 3 Fixed-fixed Bean Problem 
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and E = modulus of elasticity, A = area, I * moment of inertia. 
The authors also described the required partial derivatives such 
as the a., tensors, when variability of the material property E 

• j * 

and cross-section A or second-moment of the area I are to be 

accounted for. The "form" of the matrices remain the same. In 

the case of the calculation, the ij-th element can be 

directly obtained as i Simplified expression for 

^ can be written with a * ^ cos 2 0, b = *jr sin0cos0, d * 

6A u 

E sin 2 0 and c * e * f * g = 0, in (6.1). Similar calculations 
are possible for with a = sin 2 0, b = - sin0cos0, 

c = - ^ sin6, d « cos 2 6 , e « £§ cos0, f * fl “ 

L L L 

in (6.1) 

The aforementioned formulation does not permit the 
possibility of random description of the nodal point 


cos 2 6 , e * cos© , f * g = 
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coordinates 


In order to allow the end locations (x^,y^) and 
(xj,yj) to assume a spatial variation character, one needs to use 
in (6.2) 


e 


tan 


-1 *1 ' y i 

x . - x . 
3 i 


(6.3a) 


and 


X = ✓ , 2 , ,2 

(x i- x j + ‘yriy 


(6.36) 


This makes the 


a k 


(s) 


dx . 

l 


type formulas much more cumbersome than 


those which appear in (6.2). 

The authors describe the force vector F^ e ^ due to a 

uniformly distributed load W, which could be calculated by using 

the finite element shape functions. This entailed the 

■ (e) 


calculation of 


dF 


dW 


quantity as: 


dF 


(e ) 


dW 


a cos$ + bsin$ 

a sin$ - bcos$ - c cos4> 

a cos<t + bsin4i 

a sinij) - bcos<|> - c cos4> 


(6.4) 


L sin9 w L cos© L . 

where a = ^ — , b = ^ — , c * yy, and ♦ is the inclination of 

load W as depicted in Fig. 6.1. 

An important step in the papers is to describe the mean and 

dispersion of the nodal load due to spatial variability of the 

load distribution W(x). One utilizes the definition of the nodal 
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load : 


L 

F . - J W(x)N.(x)dx (6.5) 

■ 1 0 1 

where N^x) is the i-th shape function along x. The expected 
value operator E was applied to the above equation leading to 

L 

E [ F . ] = / E [W(x ) ] N . (x)dx (6.6) 

1 0 1 

consequently , 

L L 

E [ F . • F . ] = / / E[W(x.)W(x 2 )]N i (x 1 )N. (x 2 )dx 1 dx 2 (6.7) 

1 J 0 0 J 

Thus the mean and autocorrelation function for nodal loads were 
defined. 

The authors focused their attention on Gaussian homogeneous 
processes. If the loading function W(x) is a Gaussian process, 
then each element F^ is normally distributed. For non-Gaussian 
processes, distributions for F^ will pose computational 
difficulties. From a reliability point of view, one may have to 
restrict the computation correct up to the second-moment terms. 

The authors have developed a completely documented computer 
code FORAFS to carry out stochastical finite element analysis of 
frame structures. With the prescription of correlation 
coefficients of the basic variables, such as member area, moment 
of inertia, etc., exceedance coefficients can be calculated 
according to the first-order reliability method by using that 
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computer program. 

Liu, Belytschko and Mani [L-6.1] considered a general finite 
element representation of a dynamic system in the form: 

MU + f(X,U,U) = F ( t ) (6.8) 

where the deterministic mass matrix M is considered with the 
deterministic load vector F(t). The elastic restoring force f 
is also dependent on the state vectors: U * displacement and 

U = velocity. The paper displays very encouraging numerical 
results when compared with Monte Carlo simulation data. 

The paper makes the drastic simplification of ignoring the 
off-diagonal terms in the covariant matrix Cov[X^,Xj]. For any 
finite element system with spatial variability considerations, 
the nonzero correlation distance (which depends upon how fast the 
correlation coefficients die out, i.e., on the bandwidth of the 
correlation matrix) is indeed a key statistical consideration. 
Diagonalization of the covariant matrix indeed simplifies the 
algebra but is unrealistic for any nontrivial stochastical 
process with spatial variability. 

The algebraic derivation presented therein can be obtained 

directly from the Nakagiri and Hisada papers [N-4.1] and [N-4.3] 

when Cov[X^,Xj] is assumed in the form 6^^ Var(X^], where 

6.. is Kronecker's delta, 
ij 

The authors do not detail the simulation technique for 
nonlinear systems. As has been pointed out in various papers 
(for example, refer to [P-1.1] . Such numerical simulation is not 
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exactly a routine procedure. Secondly, the assumption of 
Gaussian distribution for spring constants as in the paper [L- 
6.1] (when the possibility of negative stiffeners is acceptable 
as a realization) is quite questionable. The paper, at best, 
could serve as an example to test out a code under the 
aforementioned stringent restrictions of covariant matrix to be a 
diagonal one. 

In this category of papers of special interest the most 
original contribution is by Contrearas, [06.1]. The state space 
representation of dynamic response for a stochastic finite 
element system was conceived to be a finite dimensional Markov 
process. The mathematical treatment is elegant and practical 
even though rather involved. Algebraic details of the paper will 
be summarized in the review of advanced methods. The following 
key steps are provided to establish a resemblance of time 
marching scheme in finite element temporal solution to a Markov 
process (where the present state depends only on the previous one 
not on the entire past history). 

The author proposed the state vector Y to house the random 
variables X, besides the usual case of displacement U and velocity 
U. Thus 



The equation of motion of a stochastical dynamic system was 
written 


90 


Y = f (Y, t] + G[Y, t] Wr 


( 6 . 10 ) 


in which f is the discrete operator describing the dterministic 
equation of motion, G is to represent the contribution of a unit 
white noie and W is a white noise. The stochastic vector 
differential equation (6.10) corresponds to Ito's form [A-l.l]. 
The following temporally discrete form was then obtained: 

5 <t n+1 ) = ♦ 5<tn>, t n+1 , tn) 

+ r{ Y( tn) , tn) w( t .. ) 

n + l 

where the discrete operators + and T were obtained in terms of 
the finite element mass, damping and stiffness matrices. 

Finally, the unknown vector Y was estimated according to Kalman 
filtering method. 


91 



7 . Outstanding Issues and Recommended Research 

The major computational concern for a successful execution of 
a finite element code with system stochasticity could be assessed 
separately for the two separate techniques, viz. the simulation 
and the perturbation (and related) methods. 

(i) Monte Carle simulation demands the execution of a 
deterministic conventional finite element code many many times . 
Each input realization (like material properties for each 
element, boundary node coordinates, temperature distribution, 
etc.) is required to be generated by a statistical package 
independent of the finite element program. A number of paper 
mainly authored by Shinozuka and his associates, [S— 7.1] , address 
the question of simulation in design-analysis for structural 
engineering problems. The theoretical background for spatially 
uniform multi-dimensional and multivariate Gaussian processes is 
quite complete. There are some computer programs available for 
research purposes which are suitable for finite element models 
with limited number of degrees-of-f reedom. There is indeed a 
need to develop robust versions of these simulation programs. 

In the case of nonGaussian and nonstationary processes 
research work is urgently required for the successful development 
of simulation algorithms. There is hardly any documented 
statistical package on the market in order to generate Monte 
Carlo database for arbitrary distributions (which could be 
prescibed in tabular observations) particularly suitable for 
discrete structural analysis. 

There is no available dynamic finite element code which is 
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integrated with a simulation procedure. Engineers, who are 
trained in computational mechanics and simultaneously possess 
working knowledge in applied statistics should be entrusted to 
develop such a statistical front end to a finite element computer 
program. 

Computers with pipeline architecture and those aided with 
parallel processors can carry out repeated deterministic 
calculations of simulation in a faster, more accurate and much 
more economical fashion, refer to [K-7.1], The Computational 
Statistics group of the Society for Industrial and Applied 
Mathematics (SIAM) organizes meetings to present the state-of- 
the-art procedures on simulation. Attention is drawn to those 
highly theoretical and analytically oriented formulations for 
Monte Carlo technique. Useful and practical computational tools 
for discrete realizations (as demanded in a finite element input 
data stream) could be developed on the basis of the 
aforementioned modern mathematical research. 

(ii) A direct finite element modeling on a stochastic input 

database invariably necessitates a Taylor series expansion at 

certain stages of computation. Thus the perturbation principle 

is quite inherent in such formulaions. In various presentations 

computation of the deviator such as &K [* K-Kg ] , refer to 

[B-1.4], or the derivatives with respect to random variables X 

such as 5K/3X^ become essential. This will demand rewriting of 

stiffness routines in a conventional finite element computer 

code. These routines are much more lengthy, especially when 

higher derivatives such as p... * 5 2 K . ./ax, 5X_ are required. 

l j eiTi 1 j x in 
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For example , in simple truss problems with stochastic nodal point 
coordinates, the expression for the derivative of the element 
stiffness matrix with respect to nodal coordinates x if y^, 

etc., involves lengthy expressions in terms of trigonometric 
inverse functions and their derivatives. The chore to hand 
calculate such derivatives and then to develop FORTRAN 
subroutines will be extremely time consuming for complicated 
finite elements, such as shell elements. 

The reviewer himself uses a computer algebra program SMP 
(Symbolic Manipulation Program), [S-7.2], to formulate finite 
element stiffness matrices and to evaluate derivatives with 
respect to the algebraic variables in the stiffness routines. 

The Taylor expansion routine in SMP is quite handy. It is 
possible to obtain readily the algebraic expressions of the 
stochastic stiffness matrix K expanded about the mean Kq when the 
closed form expression for is prescribed. It is very 

strongly recommended that attention is paid to integrate finite 
element FORTRAN program with algebraic manipulating softwares in 
order to develop versatile computer code for Probabilistic Finite 
Element Method ( PFEM ) . It is to be recognized that initial 
experience with the computational philosophy of SMP would demand 
a substantial research effort. In the long run, the code 
development activities could be expedited for PFEM and PSAM 
programs when those powerful algebraic software tools are 
implemented . 

In selecting the computer algebra program the reviewer 
recommends SMP over another similar package called MACSYMA. From 
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mathematical point of view SMP is much more versatile and is 
particularly efficient for tensor analysis. Stiffness matrices, 
and especially their derivatives constitute higher rank tensors 
which are amenable to SMP programming. Systematic and user 
friendly manipulations of those algebraic entities are more 
suited to SMP than to MACSYMA. 

The specific requirement for PFEM computer program is the 
"post-processing calculation" to estimate statistics (like mean, 
standard deviation, skewness, etc.) of response quantities (such 
as displacements, stresses, etc.) obtained according to a finite 
element formulation. The method essentially converts a continuum 
field problem into sets of matrix equations such as the strain 
(I) - displacement (U) relation: 

E a IB] 0 (7.1) 

the force (F) - displacement - (U) equation: 

F = [K] U (7.2) 

Hence the generic problem is 

Y * [A] X (7.3) 

where the statistics of X are either prescribed or computed from 
the finite element calculations. The matrix [A] is in general 
nonlinear in x. The question is then how to estimate statistics 



of Y for a large correlated [A] . The complexity of such 
expressions like (7.3) can be readily witnessed even in a linear 
elastostatic problem where the entries in the element stiffness 
matrix is a highly nonlinear function in stochastic Poisson’s 
ratio. It is a significant computational task to obtain the 
numerical values of correlation coefficients for the elements of 
stiffnes matrices, when the probability distribution function of 
the Poisson's ratio is prescribed, is a significant numerical 
task. Attention should be focused on developing a computer 
program which could yield such statistics under general 
nonGaussian prescriptions of fundamental stochastical variables. 

The final answer sought from a PFEM computer program is the 
prediction for exceedence. Reliability based methods developed 
by Wirshing and his associates (W-7.1] adequately address that 
point. Statisticians have developed Pearson series and related 
exponential families to incorporate skewness and high order 

moments in exceedence calculations. 

The adequacy of second moment based exceedence predictions 

should be examined by carefully designed Monte Carlo 
simulation. Higher order moments could contribute significantly 
in nonlinear problems especially in the case of large displace- 
ments and large strain situations. Parametric studies are recom- 
mended for benchmark problems in order to gain confidence in 
multivariate nonGaussian and nonstationary processes. Time 
history analysis for stochastic systems with a large number of 
degrees-of-freedom is essential to test the finite element codes 
generated in PSAM, PFEM programs. 
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8. 


Non-structural Applications 

Problems of both system schostacity and random forcing func- 
tion appeared in different problems of engineering science. In 
electrical and electronic communications, the noise elimination 
aspects in signal processing deal with random processes in 
time. In fluid mechanics random turbulence is of interest. 

These problems are quite closely related to those of structural 
mechanics . 

Biostatisticians in stochastic biometrics have developed 
algorithms very closely related to the stochastical finite ele- 
ments discussed herein. Bookstein [B-8.1] used both triangular 
and quadrilateral elements in a mesh to study statistical effects 
of growth in space time continuum. Goodall [G-8.1] applied 
statistical methods in a highly coupled system of nonlinear 
partial differential equations to predict plant growth. Their 
contributions are particularly important since they discuss bases 
of computational methods to solve the discrete analog with random 
system parameters. (Their work will be summarized in the 
literature review of Advanced Method.) Shinozuka and Moss- 
Salenteijn [S-8.1] applied the Monte Carlo simulation technique 
to predict growth of long bones in mammals. 

An interesting application of stochastic analysis on a 
discrete system (not really a finite element per se) deals with 
the description of branding in trees, [A-8.1], the Markov process 
which describes the stochastic behavior of the growth continuum 
degenerates into Fibonacci number series in the discrete model. 
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9 . Conferences Related to PSAM 

International Conferences on Applications of Statistics and 
Probability to Soil and Structural Engineering (Computational 
Methods dealing with probabilistic structural analysis). The 
proceedings are available in bound volumes from the host 
institutes. The first meeting [P-9.1] was in Hong Kong, 

September 13 to 16, 1971 and the fourth (last) meeting 
[p_9.2] took place in Florence, Italy, June 13-17, 1983. 

The ICOSSAR (International Conference on Structural Safety and 
Reliability) addressed the question of Structural analysis 
according to probabilistic considerations and accomodated 
analysis of discrete systems with stochastic variables. Some 
important contributions which were presented at the third 
conference, [N-4.1] and [H-4.4] on stochastic finite elements are 
reviewed here in sections 4. In the near future the fourth 
conference at Kobe, Japan, will have a session on Stochastic 
Finite Elements. The notable authors are Vanmarcke, Mochio, 
Shinozuka, Hisada, Nakagiri, and Der Kiureghian. Some of their 
research papers are reviewed in sections 4, 5 and 6. 

The National Science Foundation sponsored a recent conference 
1984, on [C-2] Water Resources where the finite element method 
was applied to nondeterministic systems. From the viewpoint of 
non-structural applications useful. These papers are useful. 

The ASCE-EMD committee on Probabilistic Methods is arranging 
a session on Stochastical Finite Element Methods at the joint 
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ASCE and ASME summer conference at Albuquerque, New Mexico, for 
June 23-26, 1985. The notable speakers are: Der Kiureghian, 

Wen, Lawrence, Ang, Shinozuka, Grigoriu, Khater and O'Rourke. 
Since 1969, ASCE also has sponsored a series of Specialty 
Conferences on stochastic mechanics and structural reliability, 
i.e. at Purdue University 1969, Stanford University 1974, Tuscon, 
Arizona, 1979 and in 1984 at U.C. Berkeley. Proceedings are 
available as ASCE publication. The majority of the papers are on 
stochastic loading rather than on system stochasticity . 

One major aspect of organizing a successful PFEM code is to 
include the state-of-the-art research in computational 
statistics. It should be remarked that engineers, who are so 
competent in computational mechanics, hardly demonstrate 
significant contribution or appreciation for the research in 
probabilistic developments. SIAM organized two conferences with 
short courses in computational statistics, the first was in 
Boulder, Colorado, June 11, 1984 and the second in Boston, 
Massachusetts, October, 21, 1984. Many of the discrete 
probabilistic formulations in PFEM have their close analog in 
^”®i3t.ed branches. Computational tools such as interactive 
graphics, use of parallel processing, employment of artificial 
intelligence in heuristic solution, would play a significant role 
in those analytical formulations. The future review on Advanced 
Methods will address the related topics such as: Nonlinear 

Optimization in Statistical Variational Formulations, Multi- 
objective Optimization Algorithms in Discrete Computations, 
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Graphical Methods in Computational Statistics, Monte Carlo 
simulations in Supercomputers and Statistical Issues and 
Uncertainty in Artificial Intelligence for discrete stochastic 
systems. Conferences organized by the American Statistical 
Association (ASA), International Association for Scientific 
Computing (IASC) and SIAM will be reviewed. 
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LEVEL 2 PFEM FORMULATION APPLIED TO STATIC LINEAR ELASTIC SYSTEMS 


1 o Index Notation for Problem Formulation 

The following discussion casts, using index notation, the Level 2 
finite element formulation for the stochastic displacement vector, element 
strains and stresses in terms of the random load vector and the random 
variables in the stiffness matrix. 

From Reference [1] the stochastic displacement vector may be written 
as 


0 * u 0 * I [k;V - k;’»ku o i 
- (k;' ak )' i K o' sF • K ; 1aKu oS 

* (k^’ak ) 2 [k-’ 4 f - r'iKuj 

. (k;' 4 k ) 3 |k;V - k‘'iku o | 


(1) 


where U Q = deterministic displacement vector 

K" 1 = inverse of deterministic stiffness matrix 
o 

AK = random stiffness matrix measured from the deterministic state 

AF = random force vector measured from the deterministic state 

The matrix is denoted as 
o 



( 2 ) 


Reference [ 1 ] : Appendix to SwRI Monthly Technical Letter to NASA-LeRC, 

date of publication: January 30, 1985. 
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which is the deterministic flexibility matrix. The size of this square 
matrix is ixj. 

The matrix AK.^ will be evaluated by expanding the stochastic stiff- 
ness matrix in a Taylor series of powers of the random variable vector 
AX about the deterministic solution. For terms up to third order in AX 


P 
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p is the number of random variables in the stiffness matrix. This 
expression for AK^ may be written in simplified index notation as 


where 


AK. . s a. ..AX + 8 . , AX AX 
ij ij* 1 l jlm l m 


+ Y. 4 . AX. AX AX 
ljtmn 1 m n 


aR. 


a . 


_L1 


l J* s! 


X 

-O 


8. 


1 


3 2 R. 


ijlm 2 


Li 


aX aX 

1 m 


X 

-o 


(U) 


(5a) 


(5b) 
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r ijimn ' 6 


3 K Li 
aX sX aX 

a m n 


(5c) 


X 

-o 


Note that values of the terms o. , B. jlm and Y iJtmn are fixed for a 
specific problem since they are the partial derivatives of the stochastic 


stiffness matrix R evaluated at the deterministic state X^. 


In equation (1), one recurring term is the matrix product K^AK. An 


element ij of this matrix is denoted as 

K" 1 AK( i j ) = K"J AK 
o J oiq qj 

which may be written using equation (4) as 


( 6 ) 


K" 1 AK( i j ) = K~ ! fa . AX„ ♦ B AX AX m 
o J oiq'- qjl l qjim l m 

+ Y AX AX AX ) 
qjitmn i, m n' 

Higher power of K^'aK may be formed using equation (7). 
(K _ 1 AK ) 2 = [K _ 1 AK(ij)][K" 1 AK(jk)] 


(7) 


For example 


( 8 ) 


where appropriate changes must be made in the subscripts in equation (7) 
We will now turn our attention to the term [K^AF - K^AKU^J in 
equation (1). The i-th component of the vector K^AF be written using 
index notation as 


k ; V(i > • Cj if j 


(9) 


110 


Likewise, the i-bh component of the vector K~^AKU becomes 

o o 

K^AKU (i) = K"! AK .U * 
o o oiq qj oj 


= k «L(°«i. aX . * 8 „... AX AX 
oiq v qjt t qjlm l m 


( 10 ) 


♦ Y n „ m JXAX AX )U . 
qjtmn t m n' oj 

The expression [k^aF - K^AKUj was denoted in equation (17) (Ref 1) 
as the vector AU^. Its i-th component may be written using equations (9) 
and (10) as 


AU.(i) = [k” 1 AF - K" 1 AKU I 
1 o o o J 


1 K oij aF J - * "qj, 


♦ X^.^AX.AX AX )u , 
qjtmn t m n' oj 


( 11 ) 


We will call equation (11), the first order term in K" 1 because it 
only involves linear terms in this matrix. Likewise, from equation (1), 
the second order term in for the vector component Al^fi) may be 
written as 


AU (i) = -(K" 1 AK) 1 (k" 1 AF - K~ 1 AKU 1 
d o ' 1 o o 0 J 


- 1 , 


,-1 


AK( iJ)AU,j ( J ) 


“K . _ (o . . AX + S , AX AX 
oiq v qjl i qjtm i m 

>Y n ,. m „AX AX AX )AU,(1) 
qjtmn t m n' 1 J 


( 12 ) 


*Note: The subscript o always represents the deterministic state. It 

does not represent an index of summation. 


Ill 



AU ^ ( j ) is evaluated from equation (11) with the appropriate subscript 
change of j for i. 

The total stochastic displacement vector 0 is the sum of the terms 

0 = U + AU, + all- + aU- + . . . ( 13) 

o 1 2 3 

where the following recursive relationship exists 


where aU^ = 

I (K _1 aF - K _1 aKU ) 
v o o o' 


au 2 . 

-(K" 1 AK)au 1 

( 14) 

au 3 = 

-(K' 1 aK)au 2 


au. = 
1 

-( K ’ 1 aK)au i _ 1 



The next part of the discussion presents an approach as to how the 
stochastic strains and stresses may be computed from the stochastic 
displacement vector. In general the element strains are computed from the 

( 3 ) 

global displacement vector and the strain-displacement matrix. If e 
represents the strain vector in element s, then 

l {s) = S (s) 0 (15) 

where is the element strain-displacement matrix. is stochastic 

if randomness can enter in the structural geometry. We will take 

g (s) = B (S) * aB (s) (16) 

o 

(si 

where B is the deterministic strain-displacement matrix. It seems 
o 

consistent with the method of computing a K [see equation (3)] to evaluate 
by expanding B^ in a Taylor series about the deterministic state. 
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This gives 


AB 


(s) . 
ij ' 


P aS- 


(s) 


LL 


1=1 aX. 


i p 

ax + •= z 

l 2 


n a 2a( s ) 

? LIll 


X 

-o 


is 1 m= i aX aX 
e m 


AX AX 
t m 


X 

-o 


n ,3a(s) 

i P P P 3 B i i 

+ t z z z ^ — 

0 t=1 ms 1 ns i aXaX aX 

& m n 


AX AX AX 
l m n 


*o 


In terms of index notation, 


(17) 


AbJ 3) S o.,AX # + B . . AX AX + y . . AX AX AX 
ij xjl t ljtm l m ljtmn l m n 


(18) 


where the o, s and y are the partial derivatives of 8^ evaluated at 

X . Therefore, a, 8, and y are entirely deterministic. 

—o 

In a similar manner, the strains can be related to the stresses. For 
no initial stresses, the relationship between the stress and strain vector 
in element s is 


^ ( s ) _ g(s)-( s) 


(19) 


where is the elasticity matrix. Since is generally stochastic 


g(s) _ c (s) + aC (s) 
o 


( 20 ) 


(s) 

where is the deterministic elasticity matrix. As in the calculation 

of AK and AB^ 3 ^ , the AC^ S ^ matrix will be computed by expanding in a 

Taylor series about the deterministic solution. Thus 


Ac( 3) : a.. a &X 0 ♦ I., B AX AX ♦ y. la AX AX 6X ♦ 
ij ljl l ijtm i m ijtmn t m n 


( 21 ) 
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_ (s) 

where a, B and y are the partial derivatives of evaluated at and 

are entirely deterministic. Obviously a relationship exists between the 

matrices AK^, AB^ and AC*^ depending on the nature of the stochastic 

state vector X^. 

Using the equations developed in this section, expressions are now 
available for the computation of the structure's random displacements 0, 
element strains and stresses entirely in terms of the random 

load vector aF and the vector AX of the random variables in the stiffness 
matrix. Obviously, the amount of computational effort depends on how many 
terms are retained in equation (13) and to what order are powers of AX 
retained in AK [equation (3)1, AB^ [equation ( 18 ) ] and AC^ [equation 
( 21 )]. 

2.0 Truncation of Stochastic Displacement Vector in Terms of Powers 
of K' 1 

Q_ 

Equation (13) expressed the total random displacement vector as 

0 = U + AU. «- AU-. ♦ AU, + . . . (22) 

o i d 5 

where AU^ involves the i-th power of the deterministic matrix K q ^. 
Successive terms for AU i are related by the factor K~^ AK. As shown by 
equation (14) 

AUi = -(K^ 1 AK)AU ._ 1 . (23) 

From the magnitudes of K* 1 and AK encountered in engineering 
problems, it is expected that AU i will be small compared to AU^. Thus, 


114 


at this stage of the formulation, it seems reasonable to take only the 
first two terms to approximate 0, i.e., 

0 = U q ♦ AU 1 (24) 

where A^ is given by equation (11). 


3.0 Truncation of Power Series Expansions for AK, AB, and AC 

The expressions for AK , AB (s ^ and AC* S * were developed in terms of a 
Taylor series about the deterministic state. These are given up to terms 
of third order by equations (4), (18) and (21). The coefficients of the 
powers of AX must be evaluated by taking the partial derivatives of the 
stiffness, strain-displacement, and elasticity matrices. 

For the present time we will confine the analysis by retaining only 
the linear terras in AX. This is a reasonable assumption if the variances 
of the random variables in the stiffness matrix remain small. Hence only 
compution of the first, and not the higher partial derivatives of the R, 
8^ and matrices is required. Under this assumption equations (4), 

(18) and (21) reduce to 


AB 


(s) 

ij 


a. <„4X. 

i j a. i 


AC 


(s) 

ij 



t 


(25) 
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^ .0 Summary of Structural Response Equations 


Under the assumptions made in the previous two sections, the 
dependent stochastic displacement, strain, and stress vectors may be 
expressed as 


0. s 11 . ♦ K~ ! .AF , - K~! a 4 „U .AX. 
i oi oij J oiq qji oj t 


(26a) 


-(s) . + ax JO. 

1 1 oij iji l J J 


(s) . - 


(26b) 


f (s) 

l 




r(s) 


(26c) 


In terms of the random variables AF and AX , the strain and stresses 
may be finally expressed as 


-(s) . ( B l s I + “ AX.] • (U + K~! AF - K~1 a U AX ] 
l 1 oij ljt t 1 1 oj ojn n ojq qnr on r 1 


.( 3 ) . - 


,-1 


,-1 


(27) 


f‘ 3) = (i 


>(s) 

oip 


5 ipi*V 


[b ( s ; 

1 opj 


a . AX 1 
pjm m 1 


(28) 



K'! AF 
ojn n 


K" ! a U 
ojq qnr on 



Equations (26a), (27) and (28) formulate the global stochastic 
structural and element response entirely in terms of the random loads and 
variables in the stiffness matrix. Note that the displacement vector 0 is 
linear in the random variables, while the strains and stresses contain 
terms up to quadratic and cubic in AF and AX , respectively. 
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10 


This formulation directly reduces to the deterministic solution if no 
random variable exist in the loading or structural stiffness. That is 



ol 



B (S >U , 
OlJ oj 


s) = c^B^U 

i oip op j oj 


(29) 


The next section applies this formulation to the specific problem of a 
three-bar truss. 


5.0 Example Applied to a Three-Bar Truss 

To demonstrate the procedures, let us consider the three-bar truss 
shown in Figure 1. The bars can take only axial loads, and displacements 
of the loaded end remain in the X-Y plane. 

The stochastic matrix equation for the system is 


RO = F 

where 



0 

x 


1 

> 

i 

II 

o 


F = 



0 


P 


yj 


yj 


The stiffness 2x2 matrix for the system is given by 


rs 


3 5 i £ i 2- 
r cos 9. 


3 JLE. 


i=1 i 


i 


Symmetric 


E - - 1 cos 0 . sin 0 . 
i=1 I. 1 1 


3 . 2 - 
E sin* 0. 


i = 1 t. 
l 


(30) 


(3D 


(32) 
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In the general stochastic finite element problem the cross-sectional 
areas (S^, moduli! of elasticity (£.), bars' lengths (* l ) and load vector 
P, may all be considered random. However, the bars’ lengths (iL) and the 
angles ((L) (see Figure 1) must satisfy geometric compatibility conditions 
even though they are random. 

In this example, we will somewhat simplify the problem by considering 
that in the undeformed state, the truss is geometrically symmetric. Thus 


9 1 = - 0 3 = 6 


(33) 


e 2 = 0 


The lengths of the bars may be random. However, the lengths of bar 1 
and 3 are related to the length of bar 2 by the compatibility relationship 


l 1 s *3 * 


COS 8 


(3*0 


The length of bar 2 (denoted now as i ) can also be expressed in 
terms of the fixed distance d and the angle 0 as 

i = d cot 6 (35) 

The stiffness matrix may now be expressed as a function of the seven 
random variables 1L, E^, and t 

.2 (s,e, - a e ] id 

< 5 1 E 1 * *3*3 > (i 2 , d 2)3/2 (J2 . d 2,3/2 


rs 




Symmetric 


(i 2 . d 2 ) 3/2 


( 36 ) 
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We will now proceed to calculate the AK stiffness matrix by the Taylor 
series expansion about the deterministic state. However, at this stage 
since we will retain only linear terms, only the first partial derivatives 
need be computed. Under this assumption 


AK = a AX 
rs rsl 1 

The terms a are evaluated from the partial derivatives of ft 

rsi rs 

partials evaluated at the determinstic solution are, after defining 

I 


(37) 


These 


(t 2 * 
1 0 

a 2 ) 1 




3ft n 
3 A ^ 

E o1 l o 

' I 3 
0 


Si. 

3E 1 * 

A ,l 2 
ol o 

I 3 

0 

ag n 

3A 2 

E o 2 
' l o 


3g 11 _ 
3E 2 = 

!o2 

t 

o 

3ft n 
3 A^ 

E 7 l 2 

o3 o 
0 


3ft n _ 
3E^ ' 

A -l 2 
o3 o 

I 3 

o 

3ft n 

(A ,E , ♦ A .E - 
. ol ol o3 o3 

)(2i 0 d 2 - 

> A o2 E „2 

3£ 


i 5 

0 


i 2 

o 

sE 12 

.Si. 

E .4 d 
ol o 

Si . 

S. # o.'o d 

3A 1 

3A 1 

t 3 

0 

3E i " 

aE t ' I 3 

o 

3ft ]2 

3A 2 

.Si. 

3A 2 

0 

3g 12 . 
3E 2 S 

«r • ° 

>R 12 

.Si. 

E ,t d 
o3 o 

3g 12 . 

3g 21 A o3 a o d 

3 A^ 

3A^ 

I 3 

0 

3E^ = 

SE 3 " I 3 

o 

3ft 12 

. Si . 

(HI,- 
*■ ol ol 

# 03 E c3» d3 ‘ 

2l 2 d ) 
o ; 

3 l 

at 


T 5 



(38) 


( 39 ) 
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s *33 

E ,d 2 
ol 

si! 33 # o/ 

3A 1 

I 3 

0 

3E 1 " I 3 
0 

•S? 

3A 2 

= 0 

= o 

3E 2 

’*33 

E -d 2 
o3 

3 ^33 A o3 d 

3A^ 

' i 3 
0 

JE 3 - j3 

o 

’*33 

3 (\,, E o, 

♦ A ,E .)d 2 l 
o3 o3 o 

31 


7* 


o 


Let us define the vector of random variables, in the stiffness matrix 
as measured from the deterministic solution, 
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Specific values of a 


now become 


rs£ 


aR 

rs 


sR 

rs 

3R 

rs 

rsl 3A 1 

a rs4 

3E 1 

a rs7 " £ 

aR 

rs 


sR 

rs 


a rs2 ' 3A 2 

a rs5 = 

3E 2 


aR 

a LI 


aR 

rs 


rs3 3A^ 

°rs6 

3E 3 



(42) 


where are evaluated from equations (38-40). 

In a linear elastic finite element formulation, the first step in 
evaluating axial strains is to compute the component of displacement 
projected on the original positions of the bars and dividing this 
displacement by the original length. Thus 


Beam 1 


Beam 2 e (2) 

Beam 3 e ^ 


cos e 0 > sin e 0 
25 1 

*1 

2 ^ 
cos 8 0^ + sin e cos 0 Cl 

i 

i 0 + d 0 
25 1 

(t 2 + d 2 ) 

0 

x 

i 

cos 0 0 - sin e 0 

25 1 




(43) 


(44) 


(45) 
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The strain displacement matrix B relates the element strains to the 


nodal displacements. 


-(S) . g(s) g 


(46) 


In this case 


i(s) 


-2 .2 
2, ♦ a 


?2 H 2 

i + d 


- d 


i 2 


-2 .2 
t d 


(47) 


5 (s) 


(48) 


In previous discussions the stochastic B was separated into 
deterministic and stochastic parts 

g(s) _ g(s) + 
o 

The deterministic part B^ 3 ^ can be evaluated by substituting i Q for 
i in equation (47). AB^ can be computed by expanding B^ S * in a Taylor 
series about I . To the first order terms in Al(AX 7 ] this gives 


AB 


(s) . 


(d 2 - l 2 )AX 7 




21 dAX- 
O 7 

-4 

1 

o 


AX. 


(d 2 - ! 2 )aX ? 

I 4 

o 


21 dAX, 
o I 


(49) 
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( s ) 

In index notation AEr ^ was expressed for element s as 

= a. . AX 
lj ljt l 


( 50 ) 


a. can be evaluated from expression (49). Thus a. =0 for all 
lj*- ljl 


ijl except 


a 117 = a 3 1 7 


^ - *o) 

1“ 

o 


21 d 

~ ~ o_ 

127 ' ' 317 ‘ r4 


(51) 


“217 = - [2 
o 

To complete the calculation of the dependent random variables, the 
expressions for the stochastic stresses must be derived. In general, 
stresses are calculated from 


- ( s) s c (s> e (s) 


(52) 


( s ) 

where C is the elasticity matrix. For the three-bar truss problem, the 
matrix only involves the modulus of elasticity. Thus, following the 
usual notation 



- 




— 


- 


E o1 

0 

0 


AE , 

0 

0 

C = C + AC = 
0 

0 

E o2 

0 


0 

4E 2 

0 


0 

0 

E 

0 3 -1 


0 

0 

AE^ 


The Taylor series expansion about the deterministic state only yields 
terms to first order since the elasticity matrix is linear in modulus. 
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In index notation AC is 


AC. 


(s) 


where AX^ = AE^ 

ax 5 = ae 2 

ax 6 = ae 3 


ij 


5. . AX 
ljt l 


This gives for a. 


ijt 


5^^ = 0 for all iji except 


(54) 


(55) 


°114 = °225 = °336 = 1 

This completes the formulation for all quantities used in equations 
(26a), (27) and (28) for calculation of the stochastic displacements, 
strains, and stresses. The next issue concerns what probabilistic methods 
can be employed to evaluate the probabilistic response of the system. 

In closing, a few comments on the question of geometric compatibility 
are in order. In this problem, the general stiffness matrix [equation 
(32)] is given in terms of the angles cL and lengths (L . The angles iL 
were eliminated in the formulation through compatibility conditions, and 
the final results only contained the random variable At(AX^). An 
alternate approach is to retain the random variables A0^ and Ai^ in the 
formulation of 0, and T^in equations (26a), (27) and (28). The 

compatibility condition could then be imposed through the correlation 
coefficients between the random variables for the probabilistic 
evaluation. 
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Three-Bar Truss 


FIGURE 1 
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where 



The mean or expected value of CL is 
E (°i) * E ( U oi * 

* E(U oi ) * Eta'jlij) (59) 

* U . ♦ a* E(a*J 

01 lj ' y 


126 


( 60 ) 


The variance of 0. is 

i 


V(0.) = E{[0 t - E(0.)] 2 ) 


= E(0. 2 ) - |E{0 1 )] ; 


which for equation (58) reduces to 

n * ,2 


V(0 ] = Z (a V(AY ) 
j=1 J J 


n n 


* » 


* 1 E “ n ® ik Cov(AY AY ) 

j = 1 k=1 lJ 1K J K 

jik 


( 61 ) 


The covariance term Cov(aY AY ) , denoted as o. 0 , is defined as 

J * J * 

CovtAYjAYj = E[ {AYj - E(AY J )}{aY k - E(AY k )}| 


= E(AY J AY k ) - E(aY,)E(aYJ 


( 62 ) 


J' 


If the random variables in 6 Y are independent, then they are 
uncorrelated and 


Cov^Yj, AYj = o Jk = 0 (63) 

In this case, the variance of 0. becomes simply 

V(0 ) = Z (a*) 2 V(AY ) (6U) 

j:1 1J J 

The correlation coefficient will not be zero in the general case and 
must be evaluated. Figure 2 illustrates how correlation could occur 
within each set of random vectors aF and AX and even between AF and AX . 
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7.0 Higher Moments and the Distribution of 0. for Independent Random 
Variables — — — 

The special case of independent, i.e., uncorrelated random variables 

permits a relatively simple calculation for the higher-order moments of 

the stochastic variable CL. Hines and Montgomery [2]* show that if 

M^ (t) is the moment generating function of AY^, then the moment 

generating function of CL for equation (58) is 

U .t 

M (t) = e 01 

0. 

i 

Recall that the moment generating function for the random variable Z is 
defined as 


l H 4Y^ a i! t )' M 4Y 2 ^ a i2 t ) 


H 4Y '“In 1 ) I 
n 


(65) 


M z (t) = E(e tZ ) 


( 66 ) 


and has the property 


d r M z (t) 


dt 



(67) 


Thus 


E(Z r ) 


d r M z (t) 

dt r 


t=0 


( 68 ) 


From knowledge of the type distributions of the independent random 
variables in AY, we know their moment generating functions. From equation 
(65), the moment generating function for (L can be constructed, and the 
r-th moment is 


m 


= E(0J) = 2— M (t) 
1 dt r 0. 

L 


( 69 ) 


1 1=0 


[2]* Probability of Statistics in Engineering Management and Science , 
Wiley, 1980. 
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The probability density function (pdf) 
moment by 


P (u.) 

0. 1 

1 


is related to the r-th 



u i p 0 

i 


( 70 ) 


In theory, at least, from a knowledge of E(O^), we can construct the 
probability density function p^ (u.) . This would allow computation of 

the cumulative distribution function (cdf) P (u.) , where the probability 

— 0. 1 
that 0. < lk is defined as 1 


P(0 s u } I P (u ) 

0 1 

1 


(71) 


U, 


J* P„ (ujdu 

. II 1 1 


8,0 £ xP ected Values, Variances a nd Covariances of the Stochastic Strain 
Vector 


From equation (27) in Section 4.0, the stochastic strain vector 
expressed in the form 


was 


J ‘ S> = 1 -J * VM'Kj * K ;>n - K oJq°qnr U on aX r1 
Expanding the expression gives 


(72) 


e i S> = B iiJ U oj * B oi j K ojn AF n ’ B oij K ojq a qnr U on AX r 


* a U*V X « * °iJt K ojn flX i ap n - ®iJt X ojq°qnr d on iX t AX r 


.-1 


-1 


(73) 
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If the product terms in the random variables are neglected in 


equation (73) , i.e. , 


AX AF * 0 
i n 

and 


( 74 ) 


AX AX * 0 
2. r 


** ( s ) 

then the strain is reduced to a linear form. Using the vector A X 

defined by equation (57), the expression (73) for strain may be written as 






(75) 


where 


.(s) . 


01 


8 <S >U , 
Olj oj 



B (s ; 

Oil o2j 


-B^f ^ K " 1 a . U 
oil olq qnj on 


+ a i2j U oi 


(76) 


The stochastic strains in equation (75) are of the same form as 
equation (58) for stochastic displacements 0^. Therefore, the development 

.( 3 ) 

of the expected values, moments, and distributions of e,. follow along 
the lines given in Sections 6.0 and 7.0. 
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9.0 Expected Values. Variances and Covariances of the Stochastic Stress 


Vector 


The same linearization strategy will be followed with the stochastic 
stress vector. From equation (28) the stress vector was expressed as 


f(s) _ [c^ s ^ + 5. ax ]-[b ( s) ♦ a ax I - 

i 1 oip ip* a 1 1 op j pjm n 1 


(77) 


[U . ♦ K“! AF - K . « 
oj ojn n ojq qnr 


-1 


U AX 1 
on r 1 


Expanding equation (77) and retaining only linear terms gives 

t| s) = C (s) B (s) .U + C^B^K' 1 AF 
1 oip opj oj oip opj ojn n 


c(s) B (s) K -1 a (j ax 
oip opj ojq qnr on r 


(78) 


* 5 ip. B op]V X . * C o!&j.V*» 


s(s) 


This expression for T. may be written in the simplified form 


• f o V * “u 4i j 


(79) 


where 


AY = 


AF 

AX 


( 80 ) 
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U> = c (») B <*)n , 

01 oip opj oj 


C . B -K"] . 
oip opt otj 



0 


0 


- C^B K -1 o .U 
oip opt otq qnj on 

♦ C (s) a U 
oip ptj ol 

+ a. ,B^ S ^U 

lpj opt ol 


( 81 ) 


The linear form of equation (79) allows the same development for 

-(s) "f ~(s) 

expected values, moments, and distributions of Tj^ as for IK and in 

Sections 6.0 thru 8.0. 


10.0 Comments on Implications of Linear Formulation 

In the previous sections the dependent stochastic variables in 
displacement, strain, and stress were cast in the linear format 


0. = 0 . + a* AY, 

i oi ij j 


; (s) 

= e (3) 

+ a* AY. 

l 

Ol 

lj J 

;j(s) 

CO • 
if— 
ii 

.« 

+ a . AY 

i 

Ol 

lj J 


(82) 


First, such a format has several nice properties in the evaluation of 
the probabilistic response. If the elements of the random variable vector 
AY are each normal and independent , then the random variables vectors 
0( s ) and are normal. Furthermore, if the elements of AYj are 

independent , then 0i s \ and approach normal distributions as 

the size of the AY vector, i.e., dimension j, becomes large. This is 
irrespective of the type of distributions in the AY vector. 
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TEMPERATURE LOAD 



o Correlation Between Material Properties 

as a Function of Location 

o Correlation Between Applied Forces 

o Correlation Between Material Properties and 

Applied Loads, e.g., Temperature Dependent 
Properties 


FIGURE 2. EXAMPLES OF CORRELATION BETWEEN RANDOM VARIABLES 
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Dr. Y.-T. Wu 
Dr. O.H. Burnside 
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September 1985 
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General Description 


3. 


4. 


,^, Th1 V\!? n d ? scr1bes te st cases designed to validate the first-year 
edition of the PSAM PFEM computer code which has the following capabilities: 

1. The code combines the conventional finite element method using the newly 
developed perturbation scheme for generating analytical performance 
function, with the fast probability Integration (FPI) algorithm. The 
performance functions will consider all Important design-related 
variables such as the displacements, strains, stresses, natural 
frequencies, limit states, etc. 

2. The code employs complex elements (beams, plane stress/straln, 
SMllTfo™it?o“ eS ' 4nd 5he " s)l ,,near e’Mtlc Mt.n»l behavior, and 

«L! tat1C a. prob1 f ms * the 1oad1n 9 vector can be treated as a correlated 

oTfriSJ^tocKW E r0bl r’ S * the stead *- state ra " d om responses 
solved' near stochastic structur e to stationary random loading can be 

°I code » for probabilistic design purposes. Include the 

funct1ons t ?rdn S nf y th UnCt1 2 nS ^ Pdf * and the cumu1at1ve distribution 
functions (cdf) of the performance functions, and the probability of 

exceedence or the reliability of the response variables. 

the PFEM codl° Wi Mo^ a nf S th ere Carefully selected to test the many features of 
tne pfem code. Most of the cases are considered well-defined, with snecifled 

numerical design values and distributional Information. However, there are 

giJen CaS Those e are°tha °h 9 °° d appr0 * imat1n 9 performance function are 
? * .f se ar f bbe cases where some trial and error procedures must be done 

to generate meaningful Input data. A case for testing probabilistic shell 
response has not been finalized yet due to the difficulty of finding f "good" 

cant11evered f nlxfp°!r f ? r ^ heck . 1n 9 Purposes. A test case involving I twined 
mJuh2 e 4 er ? d .4 P ] a ! e (f elat1n 9 closely to the turbine blade In the PSAM Drolecti 
ay be Included In the validation, depending on the availability of the 
accurate theoretical solution of the natural frequencies. 

of theVlTrnd* P ^?°! e of this test plan 1s t0 define th ® scope 

of the PFEM code validation. All test cases. In their final version will be 

analysis methods SO'?'?* S0 that the advanced reliability 

analysis methods, as well as Monte Carlo simulations, can be aoolled for the 

cX 1??« '?,?*"«">• fa " ^ thi n the^above-def 1 ned 

probabilistic deX “ C ° Ver *" "" port4nt «««* of the 


Test Cases 


Case l * Ca ntilever Beam (Solution is given in Section 4 Appendix) 

P M«| h |\ Can o 1 !f Ver beam m ° de1ed F1 gure 1 Is subjected to static loadings, 
i(1 1,5). P i s are correlated with the correlation coefficient defined as: 
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( 1 ) 


AX . . 

Op P = exp l-C(-T-^) ] 

p r j L 

where C Is a constant, ax,, is the distance from element i to element j. and L 
is the length of the beam. J J ’ L 

The tip displacement is: 


5 2P.i. c 

■ i [ (3L-a 

1=1 Ewt J 


1 - 


^ i a i L 

+ -v- 1 


( 2 ) 


where t t is the distance from the support point A to the point where P 1 

applies. The "fixed-end" is not exactly fixed, and Is modeled using a torsion 
spring with spring constant K. y 5 on 


The following distributional data is assumed: 


~ Normal (20, .1) lb 

E ~ Lognormal (10 7 , .03) psl 

L ~ Lognormal (20, .05) in. 
t - Lognormal (0.98, .05) In. 
w " Lognormal (1.0, .05) In. 

K ~ Lognormal (10 5 , .05) In-lb/rad 

o y ' Wei bull (10 4 , .10) psl 
v =0.3 

C = 1 


where X Normal (u x , C x ) means the variable X has a normal distribution 
with u x = mean value and C x » coefficient of variation. If the distribution 
is a lognormal, then u x « median. 

Another performance function considered Is the maximum stress at point A, 
denoted as S, 

s . 6lP '‘) 

j * S5 /1\ 
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The correlated loading vector, P, can be transformed to an uncorrelated 
vector, g, using 


e » A T P 


( 4 ) 


Where A Is an orthogonal matrix with column vectors equal to the eigenvectors 
of the covariance matrix, C p ; 




where 


a 1j s °1j°1 a j 


( 5 ) 


( 6 ) 


Using the Inverse transformation of Eq. (4), the two performance 
functions, namely « and S, may be written In terms of p , for example. 


a 


2 ( 3 ^ - 






Ewt' 


t-iV 


-1 


( 7 ) 


where {•} Is a column vector. The mean and the standard deviation of p can be 
computed as: 


*P 



( 8 ) 
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-I * - (9) 

where » Is e vector of the eigenvalues of C p . The deterministic solutions, 
using the mean (or median for lognormal variables) values, are 

* 0.5 In. 


and 


* 7497 psl 


which may be used to check the deterministic solutions of the PFEM code. 

The expected output of the code Include: 

1. The pdf and the cdf of the tip displacement, a. 

2. The pdf and the cdf of the maximum stress, S. 

3. The probability that the stress S will exceed the yield strength, a 

The above results of the PFEM code win be checked hv fpi nrnnram 
Monte Carlo simulation using the exact performance functions. 

Case 2. Cantilever Plate 

IJl S ^! S ca« aSe for Recking the plate element. The problems and the 
thickness !^n !!e Uken $ « S?1 In.^ U that th ® med1an 0f the 


Case 3. Cantilever Bean (Natural Frequency) 


The primary goal of this test case Is to test the capability of the 

Slv^l^Mt^! 1 ^" 1 ??!“ I he e *9 enva1ue problem. The cantilever beam, as 
?i ve ? 1n ,J est C * se w1 ] 1 be used I but the end point A will be assumed fixed 
(K~). The performance functions to be tested are the first three bend 1 no 
frequencies which may be approximated as oenaing 



, 1*1,3 


( 10 ) 


" M * ntS (a l ■ 3 - 5? - 1 * *3 ■ «•» «d • 1. th. m«„ 

® "* Lognormal (2. 5x10” 4 , 0.05) lb-$ec 2 /1n. 4 
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“ftoCj!' ( ! 0> "“ Uril freqUe " C,eS f0r «* vibration motion 

r r arsrr r-i ^ a-*? ( r- mot,on ’ the 

S?S have^pproxlmately 

By substituting the mean or median values Into Ea HOI fh* 
£?SS2! M ^? 0 S£S:. 0r "1 “ be ° bt ‘ f "' d f or ‘check 1 ng^the‘code 

“t * 497 * 9 r «d/sec, yz plane 
* 508.1 rad/sec, xy plane 

*S£g£? ° UtPUt 0f the COde ,s the « the Pdf of for both 

Case 4 ‘ -■ t * t1n ' 1 Beam Centrifugal lo a ding and Str.«s Stiffening Effects] 

,. tnl The Seometry of the beam is given in Figure 2 The tin a *1*i 
displacement due to centrifugal loading Is: 9 ** " e t1p ax1al 


2 Pq 2 z 3 


u(x-t) » § ^ 


( 11 ) 


where the variables are defined as: 

Q * 2400 rad/sec 

o ~ Lognormal (9.0X10' 4 , 0.05) lb-sec 2 /1n 4 
t - Lognormal (3.844, 0.05) In. 

E ~ Lognormal (2.9 xlO 7 , 0.1) P$ 1 

can be expressed ^v- ^ t1P d1splacement 1s ®lso a lognormal variable which 


u(x*t) - Lognormal (6.77x10 


Stress stiffening on the bendlnq 
method ( 1 1 . * 

The performance function for 
approximated as 


’ 3 , 0.188) In. 

frequency can be Included by using Galerkl 
the first bending frequency can be 


n's 
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1 


- 7^ + 2.886a 2 

Ol 

where the thickness t.ls given as 


(12) 


t - Lognormal (0.0416, 0.05) in. 

For the finite element model, the width, w. Is given as 
w - Lognormal (1.424, 0.05) in. 

By substituting a value Into Eq. (12) and letting 

H ■ (U) 

ol 

which has a lognormal distribution, Eq. (12) becomes 

w - J H + 1.662 x 10 7 (14) 

Using Eq. (14), the cdf of w,F^(u)), can be expressed in terms of the cdf of H, 

F H (h). Because H is a lognormal variable. It can be shown that 


F > 0 ) * F H (“o - 1 * 662 x 10 > 


ln(„ o 2 -1.662xl0 7 )- Uln 

* H< VS ’ 


(15) 


where * H (-) Is the standard normal cdf. Therefore, an exact distribution 
function of w can be computed easily. 

The PFEM code will be expected to generate pdf and cdf for the tip axial 
displacement and the first bending frequency. The results will be compared 
with the exact lognormal distributions. 

Case 5. Rotating Beam - Plate Element 

Repeat Case 4 using plate element. 
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Case 6. Twisted Cantilevered Plate 


A twisted plate is shown In Figure 3 where a 
a, width b and thickness h is shown clamped at one 
edge pretwisted through an angle a. 


rectangular plate of length 
edge, with the opposite 


Although much more complicated, turbine blade Is basically a twicfan 
plate and, therefore, the model shown in Figure 3 would be a oood 
test case. However, the vibration of rotation blade even for ?2. ^ ^? °? 
simple configuration considered, have found to be verv difficult * *? ,v, y 

S^p^^^he^r^ d1sagreement°has by 

resuns, especially when the aspect ratio (a/b) Is relatively small, e.g., one 

sssH^ 

perfo^ance V func“oS n b“«cu“a?e ^FclXu' U ls e “? ntU1 th4t ‘he 
included only If a suitable perfirma^^^c^cin^'dHl^d'."' ’ ' * 

CaSe 7 ‘ — te “ tth 01ff erent Corre lated Loadlnns In Different 7nn.« 

s1mp1y h supported? de The°plate r *s JXrlSdVr’"? 1 ' 7 shown ,n F ' 5ur * 4 and 1s 

having a different corre P l«ld st«?c ^Mln,"!" ,0Ur 20n “ “ Uh Hch 

The performance function considered is the center foolnt r) 
SS^SSTJ^S!^ 1 Sl " 9 ' 4 "■* " h« the fonowlh. 


6 



r z 
m*l n*l 


sin SI* sin S21 


(nr+ 




sin 


ntir 



(16) 


For simplicity, define F 
brackets In Eq. (16), l.e 


as a factor 


representing the 


series within the 


a 



(17) 


where 


0 


Et 

12(l-v 2 ) 


(18) 
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For n loads, s is 


5 


4t 

4 


n 

i 

0 1*1 


(PF) 1 


(19) 


Now consider a total of sixty-four 
loads win be divided Into four zones, 
be denoted by P^. For example, at zone 


loads as Illustrated In Figure 4. The 
| n (®®cN zone the loading vector will 


- P 1 



i 


( 20 ) 


The ioads In each zone are correlated with the correlation coefficients 
defined In each zone as 


11 

‘ * °1,16 

16,1 

* ‘ °16, 16 


( 21 ) 


The elements in the above matrix are defined as 


°U 


exp(-K 



ay 


1J 


( 22 ) 


where ax and ay are x-dlrectlon and y-dlrectlon distances, respectively. 
Using Eq. (22), a covariance matrix, C f , can be established to subsequently 
generate a transformation matrix. A,, such that 


Bi 


A T P 
-1 -1 


(23) 


where P^ Is the uncorrelated loading vector at zone 1. 

Using Eq. (23), the middle displacement can be formulated as 


a 


4,2 « 

7S A 




(24) 


which will be used as a comparison basis. 
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Because of the large amount of random variables Involved the 
of he material properties, as well as the loads and co^eUiioS col?n“en? 
in the four zones, have not been completed yet. However the list nf t-h« 8 tS 
random variables and the deterministic parameters are given as follows: 


Random Variables: E, t, v , i, Pj to P g4 
Deterministic Variables: Kj to < 4 


The goal of the analysis is to construct the distribution of «. 
Case 8. Shell Element 


Due to 
reliability 
determined. 


anaiv«i« 1C th< l8S ° f ^ e1e ^ t1n 9 4 proper performance function for 
analysis, this case (employing shell elements) has not been 


Case 9. Random Vibration f 3 1 


force T Fm « sho!n il %? P I y ^Ported beam under concentrated random 

reoresent1na S banH W iim?toH 9U wf> 5 * ^ (t) 1s a random funct1on of time, 
representing band-limited white noise with cutoff frequency w : 

c 


E[F(t)| . 0 

fS_, |w|<u 

s ,.( u .) * < ° c 

F /0, otherwise 


where E[.J denotes the mean value, S F Is the spectral density. 

as: The displacement Is a stationary random process which can be formulated 


6(x,t) * oA t Vj-V(x) ; h.( T )dt /^(c.t-T)* (Ode 
J*I J 0 J 

q(x.t) - r Qj(t)*.(x) 
j-1 J J 

-2 1 

* v j / q(x,t)Mx)dx, V 2 , » /^(xjdx 

o J J 0 J 


(25) 


(26) 


(27) 
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where hj(t) is the impulse response function associated with the jth 
mode, j (x ) is the jth mode shape function: 

3 sin ^ j-l,« (28) 


also 


o * mass density per unit area 
A * cross-sectional area * wt 


Because there is no closed form solution, the distribution of $(x t) is 
extremely difficult to construct. But $(x,t) is stationary; Its mean and 
mean-square values (in time space only) are not a function of time. In this 
test case, 

E(6(x,t)| * 0 (29) 


E[fi 2 (x,t)| 


1 i M x Hu(x)iMa)*. (a) 
j-1 k-1 J * J k 


-2 - 2 . 
"j u k 1 


Jk 


where 


(30) 


I 


jk 





t 2 2 \ 

(u j - “k> 


+ 2fl 2 ( U)j+ u> 2 ) 


(31) 


k . «•>? + - 2w,(w? - & l /A\ 

•(“ >.«*;«,)-—* — in — ^ — — ’ 

2 ("y - fi*/*) + uj + 2 u c (u* - 


1/2 

rr 


+ 20 ( tan 


- 1 


", - U - gyi)'-* 


+ tan 


-i 


0/2 

u t + (“/ ~ ^/4)' /a 


fi /2 


(32) 
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8 = 2c j W j 3 C/oA 


(33) 


where c is the damping coefficient and c. Is the damping factor of the ith 
mode* J J 


.. . joal fel lability analysis Is to construct, at x*a, the 

distribution of the mean-square value defined In Eq. (30). (This distribution 
may us ® d t0 establish the distribution of «, considering all random 


The cutoff frequency » c will be chosen such that 
"14 < w c 

where « 14 Is the mean value of the 14th natural frequency. Following is the 
list of the random variables and the deterministic parameters: 


Random Variables: E, w, t, i, c, 0 

Deterministic Value: a, S ft , u 

0 w Q 

SKCir TOdUn v, ' ues - the ,o ' iM,n9 

H - j - 0.01 
O c - w c i/(E/o) 1/2 * 2, 

B » ee/(E/ 0 ) 1/2 » 0.02 
n * | - 0.3 
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FIGURE 1. Model Definition of a Cantilever Beam 
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FIGURE 2. Rotating Beam Geometry 
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Introduction 


A preliminary plan for validation of the first year PFEM code was 
described earlier in the PSAM Monthly Progress Report No. FY '85-12 
(Attachment 4). Presented herein is the solution to the validation test 
Case 1. 

Case 1 


Problem : The cantilever beam modeled in Figure 1 is subjected to static 

loadings, Pj(i*l,5). P^s are correlated with the correlation coefficients 
defined as: 


°P 




( 1 ) 


where ax i . is the distance from element i to element j, and L is the length of 
the beam. J The goal of the analysis Is to determine: 

1. The cdf (cumulative distribution function) and the pdf (probability 
density function) of the tip displacement, «. 

2. The cdf and the pdf of the maximum stress, S, at point A. 

3. The probability that the stress S will exceed the yield 
strength, o y . 

Solution : The covariance matrix of the correlated static loadings is 



SYMMETRIC 1 


( 2 ) 


where o p = 2 is the standard deviation of P i 
The eigenvalues of C are: 


x 


14.98 
0.4385 
s 0.5941 
I 2.973 
1.059 


( 3 ) 


152 


•and the normalized matrix of eigenvectors of C p is 


0.4148 

-0.2083 

-0.3941 

-0.5871 

0.5334 

0.4621 

0.5107 

0.5871 

-0.3941 

-0.1599 

0.4782 

-0.6256 

0.0 

0.0 

-0.6163 

0.4621 

0.5107 

-0.5871 

0.3941 

-0.1599 

0.4148 

-0.2083 

0.3941 

0.5871 

0.5334 


Using Eq. 4, an uncorrelated vector is generated as 

P * A T P (5) 


because P is a normal vector, p is also a normal vector. The mean values and 
the standard deviations of p can be computed from Eq. 5 as follows. 


u 


P 




( 6 ) 


a 


P 


3.870 

0.6622 

0.7707 

1.711 

1.029 


(7) 


In Figure 1, the tip displacement is 


6 * 


Noting that 


Eq. 8 can be written as 


5 2P,., 




r I— H- (3L - l,) 

1*1 Ewt J 1 K 


is related to L by 


( 8 ) 


(9) 
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6 


(10) 


I J 

= * [0.224P. + 0.832P- + 1.728P, + 2.816 P. + 4P C ] 

Fui-h- 3 i t 3 4 5 1 


l‘ 

+ ^-.(0.2Pj + 0.4P 2 + 0.6P 3 + 0.8 P 4 + Pg) 

in which each P i can be transformed to p, using the inverse of Eq. 5 . As an 
example, Pj^ can be expressed as 

P i = 0.4148 pj - 0.2083 p 2 - 0.3941 p 3 

- 0.5871 p 4 + 0.5334 p g dD 

By substituting Eq. 11, etc., into Eq. 10, the displacement becomes a function 
of pj. 


l 3 

5 = “ 3 (4.264 p. - 0.09786 p ? + 0.3233 p^ + 2.999 p. + 0.6045 p c ) 
EWt c J ** * 


( 12 ) 


+ (1.339 pj - 0.01248 p 2 + 0.08044 p 3 + 0.6273 P 4 + 0.07842 pgj 

where the ten random variables: L, E, w, t, K and p i (1=1.5) are independent. 

The second performance function considered is the maximum stress at the 
root section. 


5 

r 

1*1 


6 i Pi 


Wl 


(13) 


* (1.2P. + 2.4P ? + 3.6 P, + 4.8 P. + 5.6 P C 1 

wt^ 1 c 3 4 5 


Using Eq. 5, Eq 13 becomes 
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( 14 ) 


S = -k (7.87 p, + 0.00856 p~ + 0.3248 p 3 
wt c J 

+ 3.529 p 4 + 0.2567 p g ] 


where S is a function of eight independent random variables. 

Note that further simplifications of Eq. 12 and Eq. 14 can be done. For 
example, Eq. 14 can be reduced to 

S = UV (15) 

where U = L/wt^ is a lognormal variable and V is a normal variable. However, 
reliability analyses using Eq. 12 and Eq. 14 provide more Information about 
all the design variables involved (e.g., the “design point" provide the 
sensitivities of the variables). Therefore, Eq. 12 and Eq. 14 are considered 
better for checking purposes, it should also be mentioned that in the PFEM 
code, a and S will be approximated by polynomial equations involving all the 
Independent random variables. The accuracies of the approximating equations 
may be checked using Eq. 12 and Eq. 14. 

The third performance function can be constructed as 

g = - s (16) 

in which a is the yield strength and S is the stress evaluated using Eq. 

14. Thus, g is a function of nine independent random variables. 

Using Eq. 12, Eq. 14 and Eq. 16, reliability analyses were performed 
using the FPI program as well as a Monte Carlo program. To check the Monte 
Carlo program, a sample size of 100,000 was used to evaluate the statistics 
of 6. The result is shown in Table 1 where the data of the ten random 
variables is listed. The median of the tip displacement Is, from Table 1, 

stimulation) = 0.40321 in. 

By substituting the medians of the random variables Into Eq. 12, the 
displacement is computed as 

5 = 0.40319 in. 

which agrees with the simulation result. 

The cdf of a, computed at thirteen values of a, is listed In Table 2, and 
is plotted on a normal probability paper (Figure 2). It shows that the 
results of the FPI analysis and the Monte Carlo simulation are close. 

The pdf of a, which is the derviative of the cdf of a, must be evaluated 
numerically using the cdf values. Therefore, for validation purposes, it is 
more direct to compare the cdf's than to compare the pdf's. However, a pdf 
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plot may be useful for making engineering judgements and decisions, therefore, 
the pdf should also be computed. For a convenient presentation of both the 
cdf and pdf, a method is suggested in the following for generating analytical 
expressions. 

Employing the one-to-one mapping: 

F(x) * *(u) (17) 


where F(x) Is the cdf of X (such as «, S) and *(u) is the standard normal cdf 
in which u Is the standardized normal variate. Assume that F(x) is known, u 
can be computed as 


u * « _1 (F(x) ] (18) 


By taking the derivatives of 


'<»> ■ ♦<“> & 


The next step is approximate 


n 

u » z a.x 
1-0 1 


1 


then the cdf and the pdf are 


F(x) » ♦( i a.x 1 ] 
i-t> 1 


and 


Eq. (17), the pdf of X is 


u by a polynomial equation: 


approximated by 


(19) 


( 20 ) 


(21) 


f(x) » *{ i 
1-0 

where the approximating 
Mathematical Functions. 


*(u) - .3989 


1 " 

a.x 1 ] • r ia.x 1-1 
1 1-1 1 

formula for *(•) is available (e.g., 
by Abramowltz and Stegun), and 

exp [-0.5u 2 ] 


( 22 ) 

in Handbook of 
(23) 
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The reason for establish Eq. 20, instead of generating 


n 

F ( x ) = z a.x 1 (24) 

i=0 1 

is because the functional relationship between F(x) and x is, in general, 
difficult to approximate using Eq. 24. On the other hand, the relationship 
between u and x is usually not significantly nonlinear. For example. Figure 2 
shows that u and 6 are approximately 1 inearly-related. Note that if X is 
normally distributed, then 

u “ * a n + a i x (25) 


In other words, u related to x linearly. 

To establish analytical expressions of F(a) and f(a), a table is 
established in the following 


5 u = * -1 [F(6)) 


0.22 

-3.649 

0.26 

-2.671 

0.30 

-1.818 

0.34 

-1.064 

0.38 

-0.391 

0.42 

0.213 

0.46 

0.764 

0.50 

1.270 

0.54 

1.736 

0.58 

2.164 

0.62 

2.564 

0.66 

2.940 

0.70 

3.290 


where the absolute values of the u's are actually the safety indices of the 
FPI output. By using a curve-fitting program, the following fourth-degree 
polynomial is established: 


4 


u * i 
1*0 


a 1 6 


(26) 


where 
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a Q = -12.1307 
a^ = -54.6602 
a 2 * -90.7163 
a 3 = 87.4329 
a 4 = -34.9088 

For the values of « computed, the relative errors In u estimates are 
approximately less than one percent. 

By substituting the coefficients of Eq. 26 into Eq. 22, it is now 
convenient to compute f($). Figure 3 shows the plot of f(6) using Eq. 22. 

The above procedure of presenting the result of the reliability analysis 
for the tip displacement can be applied for the maximum stress. 

A Monte Carlo simulation with sample size of 100,000 resulted in 


S * 7294 psi 


By substituting the medians of the random variables into Eq. 14, the stress is 
S = 7330 psi 

which is near the simulation result. 

Thirteen values of the cdf of S are computed and listed in Table 3. The 
result is also plotted on a normal probability paper (Figure 4). There is 
almost no difference between the simulation (sample size = 100,000) and the 
FPI results. 


The analytical expressions of F(s) and f(s) are given by Eq. 17 and Eq. 
19 where 


i 

u = r a.S 
i=0 1 

in which 

a Q * -12.1726 
a x * 2.88859 

a 2 - -0.245685 
a 3 = 0.0127482 

a 4 = -0.000278666 

a f(s) plot using these coefficients is shown in Figure 5. 


(27) 
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The probability that the stress exceeds the yield strength requires only 
one run of the FPI program. The result is 

p f = 0.0511 (FPI) 

« 0.0510 (Monte Carlo with sample size = 100,000) 

The pdf of the yield strength is also plotted in Figure 5 to compare with the 
pdf of the stress. The pdf of strength, which is Weibully-distributed, is 

f(x) = (f) (f) a " 1 expl-<|)“] 


where 


a = 12.0226 

8 = 10.4342 ksi 


159 



Table 1. Statistics of Tip Displacement 
for Test Case 1 

MONTE CARLO SOLUTION 

SAMPLE SIZE, K- 100000 

NUMBER OF RANDOM VARIABLES- N- 10 



RANDOM 

VARIABLES 


IABLE 

DISTRIBUTION 

KEAN/* MED I AN 

STD DEV/*CQV 

E 

LOGNORMAL 

0. 10CC0E+0B 

0. 30000E-01 

L 

LOGNORMAL 

0. 20C00E+02 

0. 50000E-01 

t 

LOGNORMAL 

0. 98C00E+00 

0. S0000E-01 

in 

LOGNORMAL 

0. 10000E+01 

0. 50000E-01 

K 

LOGNORMAL 

0. 10000E+06 

0. 50G00E-01 

Pi 

NORMAL 

0. 44643E+02 

0. 38700E+01 

p2 

NORMAL 

-0. 41530E+00 

0. 66220E+00 

P3 

NORMAL 

0. OOOOOE+OO 

0. 77075E+00 

p4 

NORMAL 

0. OOCOOE+OO 

0. 17109E+01 

p5 

NORMAL 

0. 2611 4E+01 

0. 10290E+01 


NOTE: MEDIAN AND CDV FOR LOGNORMAL VARIABLES ONLY. 


STATISTICS OF DISPLACEMENT: 
MEAN ■ 0. 40884E+00 

STD DEV - 0. 68587E-01 

MEDIAN ■ 0. 40321E+00 

COV « 0. 16776E*00 
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Table 2. The cdf of the Tip Displacement (Case 1) 


Displacement 

inch 

cdf 

a Monte Carlo 

FPI 

0.22 

0.000169 

0.000132 

0.26 

0.00468 

0.00378 

0.30 

0.0385 

• 0.0346 

0.34 

0.152 

0.144 

0.38 

0.362 

0.348 

0.42 

0.596 

0.585 

0.46 

0.785 

0.777 

0.50 

0.901 

0.898 

0.54 

0.9606 

0.9586 

0.58 

0.9858 

0.9848 

0.62 

0.9953 

0.9948 

0.66 

0.99842 

0.99836 

0.70 

0.99966 

0.999498 


a: Sample size * 100,000 
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Table 3. The cdf of the Maximum Stress (Case 1) 


Stress 

ksi 

cdf 

a Monte Carlo 

FPI 

4.4 

0.000590 

0.000599 

5.06 

0.00890 

0.00876 

5.72 

0.0545 

0.0548 

6.38 

0.188 

0.185 

7.04 

0.404 

0.400 

7.70 

0.636 

0.634 

8.36 

0.816 

0.814 

9.02 

0.922 

0.920 

9.68 

0.972 

0.9703 

10.34 

0.99107 

0.99023 

11.00 

0.99724 

0.99712 

11.66 

0.99912 

0.999219 

12.32 

0.99966 

0.999803 


a: Sample size = 100,000 
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FIGURE 1. Model Definition of a Cantilever Beam 
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Figure 2. The cdf of the tip displacement (Case 1) 
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Figure 3* PDF of Tip Displacement (Case 1) 
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Tip Displacement, inch 
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Figure 4. The cdf of the Maximum Stress (Case 1) 
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Figure 5. PDF's of Strength and Stress (C 
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Srength or Stress, ksi 
Stress ■+ Strength 
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Introduction 


In structural reliability analysis, techniques are well developed for 
computing the probability measures of a given performance function, g(X) , 
where X=(XpX 2 , ... X R ) an independent vector. For example, fast probability 

Integration (FP I ) methods are currently being used in the PSAM project. 
However, It Is not uncommon that the basic stochastic variables, X^'s, such as 
the geometrical parameters, the loadings, etc., are dependent and non-normal . 

If the joint probability distribution function of the X can be defined, 

the Rosenblatt transformation [I], as suggested by Hohenbichler and Rackwltz, 
may be employed to generate an Independent, normal vector [2j. Then the FPI 
methods may be used. Unfortunately, In practical applications, the underlying 
joint distribution functions are very difficult to construct based on either 
theory or experiment. Moreover, the Rosenblatt transformation involves the 
Inversion of the conditional distribution functions which are, in general, 
extremely difficult to derive or to compute numerically. This Is particularly 
true for the cases Involving large numbers of design random variables which 
are typical in the PSAM project. 

Another way of solving the dependency problem Is to use the marginal 
distributions and the correlation coefficients, which are relatively easy to 
obtain. It is well known that an orthogonal transformation can be employed to 
uncouple the dependency. If, in addition, the correlated variables are 
assumed normally distributed, then the transformed vector will be normal and 
Independent. On the other hand, however, if the correlated variables cannot 
be assumed normally distributed, the distributions of the transformed 
Independent vector are unknown and the FPI methods cannot be used. 

The problem associated with non-normal correlated vector has been 
addressed In the literature (3,4). For example, Der Kiureghian and Liu 
suggested that the bivariate distribution model due to Nataf (5| can be 
used. The transformed normal correlation coefficients between each two 
variables were found by Iteratively solving a double Integral equation. 

8ecause the calculation is tedious, a set of semi-empirical formulae for 
selected marginal distributions were developed. 

In this study, a different approach employing series expansion is 
developed for solving the transformed normal correlation coefficients. The 
method Is general and efficient, suitable for complicated distributions. 
Examples are provided to demonstrate the capabilities of the approximation 
method. Finally, the possible applications to the PSAM project are discussed. 

Problem Definition 


Given a non-normal vector X with marginal distributions, l.e., the 

cumulative distribution functions (cdf's) F (x.) (1=1, n), the covariance 
matrix may be constructed as X 1 1 
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where 


C 11 C 12 *•* C ln 
^nl **• C nn 


C 1j 3 E(X 1 X j 1 " E I x i lEIX jJ 


( 1 ) 


( 2 ) 


In which E(-) are the expected values. The correlation coefficients, p~ * , 
are defined as 1 J 


o 




°1°j 


(3) 


where and Oj are the standard deviations (std.) of and Xj, respectively. 

Using the bivariate normal distribution model, the normal distributions 
are established first by the following transformations: 

F x (* 1 ) = i=1 » n ( 4 ) 


where ♦(•) is the normal cdf and is a standard normal variate. Note that 

Eq. 4 defines a one-to-one mapping; therefore, x. may be formulated using the 
inverse transformation: 


X, . < 5 > 

The Inverse cdf's, i.e., F Y " l (-) are available in closed forms for 

*1 

distributions such as Weibull and extreme value. For many distribution 
models, the closed form solutions do not exist; and for a given u 
value, x must be solved Iteratively. 

The next step Is to find the correlation coefficients p^ between any u^ 
and Uj (l^j). To simplify the presentation of the analysis, consider 1=1 and 
j=2, and let p = p^ 2 (the "normal" correlation coefficient), p can be found 
by solving the following double Integration equation. 


p 


x l x 2 


1 

a l°2 



(x 1 -u 1 )(x 2 -u 2 )* 12 <lu 1 du 2 


( 6 ) 
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where * the mean values. 


*12 


1 

2ir >/~. 2 
l-o 


exp [ - 


- 2pu L u 2 
2(l-o 2 ) 



( 7 ) 


and Xj t X 2 are to be transformed to u 1# u 2 using Eq. 5. In general, there Is 

no closed form solution for Eq. 6, and the calculation of p requires 
Iteratively solving Eq. 6. The process is particularly cumbersome when the 
Inverse transformation (l.e., Eq. 5) also needs to be solved Iteratively. 

Nevertheless, If all the p's are obtained for the corresponding p Y Y , 

* 1*1 

which means that the covariance matrix of u vector is established, then an 

orthogonal transformation may be employed to construct an Independent normal 
vector suitable for reliability analysis. In the following discussion, an 
alternative procedure for computing p's will be developed which avoids the use 
of the double Integral. 

Obtaining the Normal Correlation Coefficients 
Using Series Expansion - A New Approach 

Consider two correlated random variables, denoted as Xi and X 2 . The 
correlation coefficients p Y Y can be computed as 

* 1*2 


EfXjXgJ - E[X 1 ]E[X 2 | 

"V? ‘ vi 


( 8 ) 


Define the transformation from X^ to u^ as 

X 1 = Tj(u.|) 1*1, 2 (9) 

and define 

u 2 ) * XjX 2 * T^(u^)T 2 (u 2 ) (10) 

Eq. 8 may be expressed as 

p X|X 2 ° 1°2 * T i (U) 
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Expand H into Taylor series about the point (uj, U 2 ) = (0,0), 


H * H 00 + 2^ H 20 U 1 + 2H 11 U 1 U 2 + H 02 U 2^ + H ’°* T ' 


( 12 ) 


where H.O.T. * Higher Order Terms 




H 


Mj 


au jau 2 


u=0 


d^Tj 


du 


d J T, 


du: 


(13) 


u=0 


Eq. 11 can be written as 

°X 1 X 2 °1 0 2 = H 00 + 2 H 20 E ^ u 1* 
- ElTj] - E(T 2 1 


2H 11 E{u 1 u 2 1 + H 02 E * U 2* + E(H * 0 * T| 

(14) 


The expected values of uju^ can be derived using the moment generating 
function. For example, it can be shown that: 


E[u 1 u 2 1 = o 
Eluf 1 = 1 
E[u|) = 0 

E[u 2 u 2 | = 1 + 2 o 2 
etc. 


(15) 


A derivation of E[uju^| up to i+j*8 has been carried out and Is listed In 
Table 1. By using Eq. 15, Eq. 14 simplifies to 



Applying the condition that o=0 when p y Y =0, Eq. 16 further reduces to 

* 1*2 


X a l°2 3 ^^11 ^ ElH.O.T. | 


(17) 


Therefore, a rough approximation, neglecting the fourth and higher-order 
(third-order coefficient is zero) terms, is simply 


o 


°l q 2 

H 11 


%X 


12 


(18) 


in which 


H 


11 


dT, 


dUi 



u*0 


(19) 


where defines the transformations, e.g., Eq. 5. 

If T i are linear, e.g., when X i are normal distributions, with 
means and standard deviations o^, then 

X 1 * *1" -j (*J ) “ + o 1 (20) 


Employing Eq. 19 gives 


H 11 = °l a 2 


( 21 ) 


and Eq. 8 then degenerates to 


p 3 



( 22 ) 


as expected. 

In general, are non-linear functions of Ui, therefore, Eq. 18 is a good 
approximation only when X 1 are not significantly non-normal. Using Table 1, a 
more complete approximation formula, up to eighth-order terms of H series was 
derived as: 
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°X l X 2 0 l°2 = o{H ll + ?( H 13 +H 3l) + 5 (H L5 +2H 33 +H 51^ + ^ H 17 +3H 35 +3H 53 +H 71 ) 1 

+ |-[H 22 + |( H 24 +H 42) + ^ H 26 +2H 44 +H 62^ 

3 . (23) 

+ f“ [H 33 + 2^ H 53 +H 35^ 

+ 24 ^ 44 ' 

This equation is believed to be adequate for the practical problems involving 
highly non-linear transformations. The inclusions of even higher-order terms 
are straight forward using Eq. 23. 

The procedure of computing p may be outlined in the following steps, 
suitable for computer programming: 

1. Select N (say N=9) points of u values, e.g., from -4 to +4, with 
increments of 1. Compute X^*T(u^) (1*1, 2) for the N values of u. 

2. For both T 1 and T 2 , find the (N-l)th order approximating polynomials 
using proper numerical schemes. 

3. Compute using the result of step 2. 

4. Solve o from Eq. 23. 

For small coefficient of variations (COV's), say COV < 0.15, typical of 
the engineering problems, the fourth-order approximation provides a relatively 
efficient way of computing p. The equation is 

0 X 1 X 2 °l , ’2 * »l H ll4 (H 31* H 13 )l * f* H 22 (24> 

The highest derivatives for X 1 and X 2 are the third orders and p may be found 
by solving a quadratic equation. 

Examples 

To demonstrate how to use the derived formula, and to examine the 
accuracy, consider the following examples. 

Example 1 - Xj^ and X 2 are both lognormal ly distributed, with equal COV's * 

C. The transformations are: 
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(25) 


inX j “ u y 


F x/ X i^ = — ) = 


or 


X* 3 6Xp[yy + U *Ow ] 

i T i 1 M 


(26) 


where Y. = mX^ 1=1,2) and 


Uy 2 3 mX 4 * in 


— fcllrt * — lr I I - 

' ST?'' 


(27) 


Jy^ ««/ m(l+C^) * ffy 


(28) 


where X^ Is the median of X^ 


From Eq. 26 


d"x 


— °y exp K + u i a v ] 

du? Y 1 Y i 1 Y 


(29) 


It follows that 


H 1j * (® Y ) 1+ '* ex P[» 1 Y 1 + »*Y 1 


(30) 


Substituting Eq. 30 Into Eq. 23 without truncating the series, and using Eq. 
27, It can be shown that 


2 2 

° X 1 X 2 C * eXp ^ po Y^ _1 


( 31 ) 
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Therefore, 


2 

tn(l+P Y y C ) 
* 1*2 

m(l+C 2 ) 


(32) 


which Is an exact solution. The more general solution for the case C^C 2 
derived using a similar approach is 


ln ^ 1+0 X 1 X 2 C l C 2^ 

/ in(l+C 2 )tn(l+C 2 ) 


(33) 


In general, closed-form solutions may be very difficult to derive or 
simply doesn't exist. In this case, a proper truncating series may be used. 
In the following, the effect of nonlinear transformation and the effect of 
truncating the series will be discussed using the two lognormal variables 
case. 

Recall that the nonlinear transformation of a lognormal variable X to a 
normal u is 


X 


Xexp[utn 



(34) 


Ignoring the constant X, the functional relationship between X and u is 
plotted in Figure 1 for three values of C, namely C=0.1, C=0.3 and C=0.5. It 
shows clearly that the C values significantly affect the nonlinearity of 
X*T(u) around u=0. 

Now consider the truncated series. Assume that Cj * 0.1 and C 2 = 0.5, 
the exact solution is: 


p * 21.222 *n(l+.05p ¥ Y ) (35) 

* 1*2 

and the approximating solution, using Eq. 23, is 
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.05618o x ^ = O [.04712+.005492+.00032+. 000012431 
+ p^[ .00111+. 0001294+. 00000754) 

(36) 

+ p J (.00001744+. 000002031 


+ p 4 (. 000000205) 


where It Is evident that the series converges quickly. As an example, 
let Py y — 0.9 , the results are as followsi 
* 1*2 

p (Second-order approx.) = 1.192 

p (Fourth-order approx.) = 0.9423 

p (Sixth-order approx.) * 0.9345 

p (Eighth-order approx.) = 0.9341 

o (Exact) = 0.9341 

It may be observed that (for p Y > o case) 

*l x 2 


° 1°2 

H U \x 2 > p exact > PXjXg (37) 


which means that the second-order approximation (using and only) 

provides an upper bound of exact p. Note that |p| (Exact)>|p Y Y I has been 

* 1*2 

proven (e.g., see Lancaster (6)). An Important application of Eq. 37 is that 
If the bounds are judged narrow enough, then it Is not necessary to try to 
obtain very accurate p value. 

Examgle2. Consider the case where one (say Xi) of the two variables is 
normally distributed, then 1 


dX l 

■ A- a O 

1 (38) 
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(39) 


d n X 


— = 0 for n 


du 


> 1 


at u^*0. 


The approximating series of Eq. 23 simplifies to 
‘»X 1 X 2 °2 “ ^7 (H 11 + ? H 13 + i H 15 + H 17 1 


%X, 2 oj 


3 

dx~ « d 

» [ dur*?-r * " l 


du: 


U2=0 


(40) 


where p Is not a function of X^. 


Assume that X ? is a lognormal variable with median of unity and COV of 
0.5. Pretend thavthe differentiation of x, with respect to u, Is difficult 
and therefore must be done numerically. Using the strategy suggested earlier, 
nine sets of solutions are obtained as follows: 


Set 

u 2 

T(Uj) 

1 

-4 

0.1511 

2 

-3 

0.2424 

3 

-2 

0.3887 

4 

-1 

0.6235 

5 

0 

1.0 

6 

1 

1.6038 

7 

2 

2.5722 

8 

3 

4.1253 

9 

4 

6.6162 


The function T(u 2 ) was plotted In Figure 1. 

The next step Is to construct an eighth-order polynomial denoted as 
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8 

x- = Z A u!J 
2 n*0 n 2 


(41) 


The required derivatives for computing 0 are 


d x. 


du. 


U2 = 0 


= A • n 
n 


(42) 


where n«l, 3, 5 and 7. 

By solving nine simultaneous linear equations, the coefficients are found 

dS • 


= 0.472353 
A 3 - 0.017595 
A 5 - 0.000191 
A 7 » 0.000013 

Using Eq. 42 and Eq. 40, the approximation solution is 


o * 1.0584 o 


X 1 X 2 


The exact solution can be derived as: 

^2 

Ox % * 1.0584 o 


y tn(l+c|) 1 2 


X 1 X 2 


(43) 


(44) 


which proves that the proposed algorithm works very well. 

Applications 

Before discussing the possible applications. It Is worthwhile to note 
that: (a) the correlations between the design variables may have significant 
effect on structural analysis (e.g., see Thoft-Chrlstensen [6]), and (b) In 
probabilistic finite element analysis, the loading as well as the geometry 
must be discretized. For small element size, correlations usually exist 
between adjacent elements. 
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Let us consider a long bar. The cross-sectional area may be treated as a 
random variable. By discretizing the bar into n elements, there are n areas, 
Ai, each of which Is a random variable. If the element lengths are relatively 
snort, then It would be unrealistic to assume that the adjacent are 

independent because the Independency suggests the sudden changes In areas. On 

the other hand, A^s cannot be assumed perfectly correlated If, in reality, 
A^'s are changing along the bar. 

A possible solution to this problem is to treat the area (along the bar) 
as a random process. If the bars are manufactured under quality control, then 
it seems reasonable to further assume that the random process Is stationary. 
Under the above assumptions, the marginal distribution as well as the 
correlations may be extracted from the measurement data. Obviously, the area 
need not be normally distributed. 

The above discussion may be extended to two dimensional problems. For 
example, the thickness of a nominally flat plate may be treated as a 
stationary random process. The correlation functions may be constructed along 
different directions. This approach Is, In fact, very similar as in defining 
a correlation field of random loading. 

The treatment of the material property, e.g., modulus of elasticity 
(again, may be non-normal) as a correlated; but not perfectly correlated, 
random field Is much more involved because it would be difficult to obtain 
experimental data for small elements. Correlation function needs to be 
assumed or established using other material properties which are related to 
the strength of material (e.g., Brlnell Hardness Number). For the PSAM 
project, the selections of the correlations need to be tailored to the 
specific problem under investigation. For example, the modulus of elasticity 
of a turbine blade may be considered perfectly correlated. However, some 
Independency may be assumed among different blades. 

Assuming that the correlation functions are defined for the elements, the 
correlated non-normal variables can be transformed to Independent normal 
variables using the procedure proposed earlier. Because the number of 
variables may be large, it is suggested that the orthogonal transformation 
should be done on a zone-by-zone basis where a zone Is defined as a region in 
which the variables are correlated; there is negligible correlation outside 
this region. Since the number of variables in each zone is relatively small, 
the computation time may be reduced significantly. 

Summary 

Given a non-normal dependent vector with marginal distributions and 
correlation coefficients, a method using series expansion was developed for 
obtaining approximating correlation coefficients of the transformed joint 
normal distributions. Then the orthogonal transformation may be Implemented 
to obtain independent normal vectors suitable for FPI or other reliability 
analysis. 

Several levels of approximations were obtained and their accuracies and 
usefulness discussed. For example, the second-order approximation (using only 
Hu), which is easy to obtain, may provide a close bound of the exact 
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solution. The series expansion has been derived up to eighth-order, which 
should be adequate for the problems encountered In the PSAM project. A simple 
numerical algorithm was suggested for computer programming. 

Finally, In discussing the applications of the developed method. It was 
suggested that the geometries, the material properties, etc., may be treated 
as non-normal dependent vectors, and that the orthogonal transformation should 
be performed on a zone-by-zone basis. 
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Table 1 


Expected Values of the Functions of Two Correlated 
Standard Normal Variables (u, v) 


FUNCTION 

ORDER 

EXPECTED VALUE 

U, V 

2 

0 

u 2 

2 

1 

uv 

2 

0 

u^.v^.u^v.uv 2 

3 

0 

c 

< 

4 

3 

u^v.uv 2 

4 

3o 

U 2 V 2 

4 

1+2 o 2 

uV.uVuv 4 . etc. 

5 

0 

u 6 ,v6 

6 

15 

5 5 

LTV, uv 3 

6 

15o 

U 4 V 2 ,U 2 V 4 

6 

3+12p 2 

U 2 V 2 

6 

9p+6p 3 

u^,v^,u®v,uv®, etc. 

7 

0 

u®,v® 

8 

105 

u^v.uv^ 

8 

105p 

U®V 2 ,U 2 V® 

8 

15+90p 2 

U 5 V 2 t U 3 V^ 

8 

45p+60p 3 

U 4 V 4 

8 

9+72p 2 +24p 4 
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Figure 1. Nonlinear Transformation of a Lognormal Variable 
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